AFRL-RQ- WP-TR-20 1 4-0051 


ENERGY-BASED  DESIGN  OF  RECONFIGURABLE 
MICRO  AIR  VEHICLE  (MAV)  FLIGHT  STRUCTURES 

James  J.  Joo  and  Gregory  W.  Reich 

Design  and  Analysis  Branch 
Aerospace  Vehicles  Division 

Richard  V.  Beblo 

University  of  Dayton  Research  Institute 


FEBRUARY  2014 
Final  Report 


Approved  for  public  release;  distribution  unlimited 


AIR  FORCE  RESEARCH  LABORATORY 
AEROSPACE  SYSTEMS  DIRECTORATE 
WRIGHT-PATTERSON  AIR  FORCE  BASE,  OH  45433-7542 
AIR  FORCE  MATERIEL  COMMAND 
UNITED  STATES  AIR  FORCE 


NOTICE  AND  SIGNATURE  PAGE 


Using  Government  drawings,  specifications,  or  other  data  included  in  this  document  for  any 
purpose  other  than  Government  procurement  does  not  in  any  way  obligate  the  U.S.  Government. 
The  fact  that  the  Government  formulated  or  supplied  the  drawings,  specifications,  or  other  data 
does  not  license  the  holder  or  any  other  person  or  corporation;  or  convey  any  rights  or 
permission  to  manufacture,  use,  or  sell  any  patented  invention  that  may  relate  to  them. 

This  report  was  cleared  for  public  release  by  the  USAF  88th  Air  Base  Wing  (88  ABW)  Public 
Affairs  Office  (PAO)  and  is  available  to  the  general  public,  including  foreign  nationals. 

Copies  may  be  obtained  from  the  Defense  Technical  Information  Center  (DTIC) 
(http://www.dtic.mil). 

AFRL-RQ-WP-TR-20 14-0051  HAS  BEEN  REVIEWED  AND  IS  APPROVED  FOR 
PUBLICATION  IN  ACCORDANCE  WITH  ASSIGNED  DISTRIBUTION  STATEMENT. 


*//Signature// 


//Signature// 


JAMES  J.  JOO,  Engineer 
Design  and  Analysis  Branch 
Aerospace  Vehicles  Division 


THOMAS  C.  CO,  Branch  Chief 
Design  and  Analysis  Branch 
Aerospace  Vehicles  Division 


//Signature// 

FRANK  WITZEMAN,  Chief 
Aerospace  Vehicles  Division 
Aerospace  Systems  Directorate 


This  report  is  published  in  the  interest  of  scientific  and  technical  information  exchange,  and  its 
publication  does  not  constitute  the  Government’s  approval  or  disapproval  of  its  ideas  or  findings. 

*Disseminated  copies  will  show  “//Signature//”  stamped  or  typed  above  the  signature  blocks. 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite 

1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YY) 

February  2014 

2.  REPORT  TYPE 

Final 

3.  DATES  COVERED  (From  ■  To) 

01  October  2010-30  September  2013 

4.  TITLE  AND  SUBTITLE 

ENERGY-BASED  DESIGN  OF  RECONFIGURABLE  MICRO  AIR  VEHICLE 
(MAV)  FLIGHT  STRUCTURES 

5a.  CONTRACT  NUMBER 

In-house 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

61102F 

6.  AUTHOR(S) 

James  J.  Joo  and  Gregory  W.  Reich  (AFRL/RQVC) 

Richard  V.  Beblo  (University  of  Dayton  Research  Institute) 

5d.  PROJECT  NUMBER 

3002 

5e.  TASK  NUMBER 

N/A 

5f.  WORK  UNIT  NUMBER 

Q0CU 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Design  and  Analysis  Branch  (AFRL/RQVC) 

Aerospace  Vehicles  Division 

Air  Force  Research  Laboratory,  Aerospace  Systems  Directorate 
Wright-Patterson  Air  Force  Base,  OH  45433-7542 

Air  Force  Materiel  Command,  United  States  Air  Force 

University  of  Dayton 
Research  Institute 

300  College  Park 

Dayton,  OH  45469 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFRL-RQ-WP-TR-20 1 4-005 1 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 

Aerospace  Systems  Directorate 

Wright-Patterson  Air  Force  Base,  OH  45433-7542 

Air  Force  Materiel  Command 

United  States  Air  Force 

10.  SPONSORING/MONITORING 
AGENCY  ACRONYM(S) 

AFRL/RQVC 

11.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER(S) 

AFRL-RQ-WP-TR-20 1 4-005 1 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 


13.  SUPPLEMENTARY  NOTES 

PA  Case  Number:  88ABW-2014-0632;  Clearance  Date:  20  Feb  2014. 

14.  ABSTRACT 

The  objective  of  the  project  is  to  understand  how  to  mechanize  multi-jointed  MAV  wings  for  perching  and/or  flapping 
applications  and  develop  an  energy-based  design  framework  for  the  solution  of  combined  multi-physics,  multi-objective 
problems. 


15.  SUBJECT  TERMS 

MAV 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT: 

SAR 

18.  NUMBER 

OF  PAGES 

102 

19a.  NAME  OF  RESPONSIBLE  PERSON  (Monitor) 

James  J.  Joo 

19b.  TELEPHONE  NUMBER  (Include  Area  Code) 

N/A 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39-18 


TABLE  OF  CONTENTS 


Content  Page 


List  of  Figures . ii 

List  of  Tables . iv 

Preface . v 

Section  1.  Introduction . 1 

1 . 1  Morphing . 1 

1 .2  Optimization . 4 

1 .3  Problem  Description . 5 

Section  2.  Methodology . 7 

2. 1  Geometry  Setup . 7 

2.2  Finite  Element  Model . 12 

2.2.1  Step  1:  Discretization . 13 

2.2.2  Element  Stiffness  Matrix . 13 

2.2.3  Global  Stiffness  Matrix . 26 

2.2.4  Displacement  Solution . 27 

2.2.5  Element  Stresses . 28 

2.3  Aerodynamics . 29 

2.3.1  Perching  Maneuvers . 29 

2.3.2  Force  Estimation . 32 

2 .4  Optimization . 35 

2.4. 1  Compliance  Objective  Using  the  SIMP  Model . 35 

2.4.2  Optimality  Criteria  Method . 40 

2 . 5  Substructure  and  Actuation  Optimization . 42 

2.5.1  Elements . 43 

2.5.2  Element  Formulation . 44 

2.5.3  Vortex  Lattice  Method  by  Ring  Vortex  Elements/Aeroelastic  Effects . 49 

2.5.4  Optimization  Background . 51 

2.5.5  Problem  F  ormulation . 57 

2.5.6  Sensitivity  Analysis . 59 

Section  3 .  Results . 61 

3 . 1  Aerodynamic  Results . 61 

3 .2  Optimized  Static  Structural  Layouts . 67 

3 . 3  Grid  Independence . 82 

3.4  Optimized  Actuation  System . 83 

Section  4.  Conclusions . 85 

4.1  Summary . 85 

4.2  Recommendations . 86 

4.3  Future  Work . 87 

List  of  Acronyms,  Abbreciations,  and  Symbols . 88 

Section  5.  Bibliography . 89 


Approved  for  public  release;  distribution  unlimited. 


LIST  OF  FIGURES 


Figure  Page 


Figure  1:  Eagle  Owl  Soaring  [2] . 1 

Figure  2:  Eagle  Owl  Flyby  [3] . 2 

Figure  3:  Eagle  Owl  Perching  [4] . 2 

Figure  4:  Comparison  of  Fixed-Wing  and  Morphing  Performances  [5] . 3 

Figure  5:  Lockheed-Martin’s  Z-Wing  [6] . 4 

Figure  6:  NextGen’s  Batwing  [6] . 4 

Figure  7:  Generalized  Topology  Design  Problem . 5 

Figure  8:  Geometry  Layout  for  Three  Section  Wing . 7 

Figure  9:  Wing  Geometry  GUI . 8 

Figure  10:  Example  Wing  Geometry . 9 

Figure  1 1 :  Birdwing  in  Forward  Swept  Configuration  (m) . 9 

Figure  12:  Birdwing  in  Zero  Sweep  Configuration  (m) . 10 

Figure  13:  Birdwing  in  Back  Swept  Configuration  (m) . 10 

Figure  14:  Birdwing  in  Dive  Configuration  (m) . 1 1 

Figure  15:  General  Quadrilateral  Element  Mapping  to  Computational  Domain  [13] . 14 

Figure  16:  Mapping  of  General  Quadrilateral  Element  to  Computational  Domain . 19 

Figure  17:  Natural  Coordinate  Unit  Vectors . 25 

Figure  18:  Computational  Domain  Indexing  Example . 27 

Figure  19:  Computational  Domain  Indexing  Example . 28 

Figure  20:  Perching  Maneuver:  Altitude  vs.  Range . 30 

Figure  21:  Perching  Maneuver:  Range  vs.  Time . 30 

Figure  22:  Perching  Maneuver:  Attitude  vs.  Range . 31 

Figure  23:  Perching  Maneuver:  Velocity  vs.  Range . 31 

Figure  24:  Main  Menu  of  Tornado  Software . 32 

Figure  25:  Example  of  Tornado  Geometry  Output  [21] . 33 

Figure  26:  Example  of  Tornado  Vortex  Panels  Output  [22] . 33 

Figure  27:  Viscous  Drag  Estimation  Curve . 34 

Figure  28:  SIMP  Penalization  of  Element  Thickness . 38 

Figure  29:  Example  “SIMP”  Output . 41 

Figure  30:  Example  “SIMPm”  Output  with  No  Checkerboard  Filter . 42 

Figure  31:  Model  Element  Designation . 43 

Figure  32:  Line  Element  Configurations . 43 

Figure  33:  Joint  Element  Configurations . 44 

Figure  34:  Typical  Joint-Beam-Joint  Model . 46 

Figure  35:  Membrane  Pre-Strain . 49 

Figure  36:  Cp  Distributions . 50 

Figure  37:  Static  Aeroelastic  Analysis . 51 

Figure  38:  Cp  for  Birdwing  Geometries  along  Perching  Trajectory . 63 

Figure  39:  Force  [N]  for  Birdwing  in  Swept  Back  Configuration  at  Point  1 . 64 

Figure  40:  Force  [N]  for  Birdwing  in  Dive  Configuration  at  Point  2 . 65 

ii 

Approved  for  public  release;  distribution  unlimited. 


Figures  (List  of  Figures  continued)  Page 


Figure  41:  Force  [N]  for  Birdwing  in  Zero  Sweep  Configuration  at  Point  3 . 66 

Figure  42:  Force  [N]  for  Birdwing  in  Forward  Swept  Configuration  at  Point  4 . 67 

Figure  43:  Membrane  Structure  for  Birdwing  at  Point  1 . 68 

Figure  44:  Bending  Structure  for  Birdwing  at  Point  1 . 69 

Figure  45:  Combined  Structure  for  Birdwing  at  Point  1 . 70 

Figure  46:  Combined  Structure  for  Birdwing  at  Point  1  without  Viscous  Drag . 71 

Figure  47:  Membrane  Structure  for  Birdwing  at  Point  2 . 72 

Figure  48:  Bending  Structure  for  Birdwing  at  Point  2 . 73 

Figure  49:  Combined  Structure  for  Birdwing  at  Point  2 . 74 

Figure  50:  Combined  Structure  for  Birdwing  at  Point  2  without  Viscous  Drag . 75 

Figure  51:  Membrane  Structure  for  Birdwing  at  Point  3 . 76 

Figure  52:  Bending  Structure  for  Birdwing  at  Point  3 . 77 

Figure  53:  Combined  Structure  for  Birdwing  at  Point  3 . 78 

Figure  54:  Combined  Structure  for  Birdwing  at  Point  3  without  Viscous  Drag . 78 

Figure  55:  Membrane  Structure  for  Birdwing  at  Point  4 . 79 

Figure  56:  Bending  Structure  for  Birdwing  at  Point  4 . 80 

Figure  57:  Combined  Structure  for  Birdwing  at  Point  4  without  Viscous  Drag . 81 

Figure  58:  Combined  Structure  for  Birdwing  at  Point  3  with  Coarse  Mesh . 82 

Figure  59:  Convergence  Data  for  Figure  58 . 83 

Figure  60:  Results  without  Membrane  Elements . 83 

Figure  61:  Results  with  Membrane  Elements . 84 


iii 

Approved  for  public  release;  distribution  unlimited. 


LIST  OF  TABLES 


Table  Page 


Table  1:  Sectioned  Geometry  for  Multiple  Configuration  Birdwing . 12 

Table  2:  Data  for  Selected  Points  along  Perching  Maneuver . 32 

Table  3:  Aerodynamic  Data  for  Birdwing  along  Perching  Trajectory . 61 


iv 

Approved  for  public  release;  distribution  unlimited. 


PREFACE 


A  topology  optimization  model  for  conceptual  wing  structure  layouts  of  morphing  micro  air  vehicles 
(MAVs)  has  been  developed  and  implemented  in  MATLAB.  Specifically,  a  6  degree-of-freedom  (DOF) 
finite  element  (FE)  model  with  a  general  quadrilateral  discretization  scheme  was  created  by  superposition 
of  a  known  simple  linear  plane  membrane  element  and  a  Kirchhoff  plate  bending  element  derived  herein. 
The  purpose  of  the  six  degree-of-freedom  model  was  to  accommodate  in-plane  and  out-of-plane 
aerodynamic  loading  combinations.  The  FE  model  was  validated  and  the  MATLAB  implementation  was 
verified  with  classical  beam  and  plate  solutions.  A  compliance  minimization  optimization  objective  was 
then  formulated  with  the  Solid  Isotropic  Material  with  Penalization  (SIMP)  method,  subject  to  the 
equilibrium  constraint  computed  by  the  FE  model,  and  solved  with  the  Optimality  Criteria  (OC)  method. 
With  the  topology  optimization  model  in  place,  four  aerodynamic  loading  scenarios  were  extracted  from 
points  along  a  feasible  MAV  perching  flight  trajectory  and  used  to  determine  wing  thickness  distributions 
for  given  planform  shapes.  The  results  suggest  conceptual  structural  layouts  in  morphing  MAVs,  but 
equally  important,  the  simple  MATLAB  implementation  of  the  model  can  be  adapted  for  a  variety  of 
objective  statements  for  MAV  morphing  wing  design. 
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Section  1.  INTRODUCTION 


1.1  Morphing 

The  ambiguous  characterization  of  an  aircraft  as  morphing  could  denote  any  one  of  sundry  possible 
modifications  to  an  aircraft’s  structure  during  flight;  but  a  more  recent  goal  in  this  venue  seeks  to 
characterize  morphing  aircraft  as  those  that  can  significantly  modify  their  wing  shape  to  adapt  to  multiple 
mission  roles,  hence  developing  high-performance  in  dissimilar  flight  regimes. 

In  2003,  the  Defense  Advanced  Research  Projects  Agency  (DARPA)  began  the  2  1/2  year  Morphing 
Aircraft  Structures  (MAS)  program  in  conjunction  with  the  Air  Force  Research  Laboratory’s  (AFRL)  Air 
Vehicles  Directorate  (now  Aerospace  Systems  Directorate).  The  MAS  program  sought  to  design  and 
build  aircraft  where,  distinct  from  all  of  the  historical  morphing  efforts  of  the  past,  morphing  would  entail 
radical,  aerodynamically  efficient  shape  change,  thus  transforming  the  mission  of  an  aircraft.  DARPA 
further  elaborates  that  “the  program  envisions  changing  wing  areas,  wing  spans,  and  other  dimensions  far 
more  radically  than  before,”  where  “radical”  describes  a  change  on  the  order  of  50  percent.  [1]  Depending 
on  the  mission  requirements,  these  efforts  will  likely  convert  military  aircraft  from  large  and  heavy, 
piloted  aircraft  to  relatively  small,  unmanned  aerial  vehicles  (UAVs).  That  nature  does  in  fact  take 
advantage  of  drastic  wing  morphing  through  disparate  flight  maneuvers  is  piquantly  and  elegantly 
epitomized  by  the  eagle  owl  in  Figure  1  to  Figure  3,  in  which  the  planform  shapes  are  outlined. 


Figure  1:  Eagle  Owl  Soaring  [2] 
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Figure  2:  Eagle  Owl  Flyby  [3] 


Figure  3:  Eagle  Owl  Perching  [4] 


Modifying  drag  is  crucial  as  lessening  drag  is  a  priority  for  both  long  range  endurance  missions  and  high 
speed  maneuvers,  and  accumulating  drag  is  paramount  for  rapid  deceleration  of  a  vehicle  to  a  landing  of 
near  zero  velocity.  But  producing  a  change  in  either  lift  or  drag  should  not  be  thought  of  as  separate 
affairs;  rather  aerodynamic  efficiency  is  influenced  each  time  the  lift-to-drag  ratio  is  changed.  The  “best” 
velocity  for  a  particular  mission  segment  often  depends  directly  on  wing  plan-form.  An  effective  graphic 
juxtaposing  conventional  fixed  geometry  and  morphing  geometry  is  that  of  the  spider  plot  shown  in 
Figure  4.  In  the  plot,  the  outermost  radius  of  the  plot  is  the  best  possible  performance  for  a  particular 
mission  segment.  On  the  legend,  the  “Firebee”  represents  the  performance  of  a  fixed-wing  aircraft,  the 
“airfoil”  represents  the  same  aircraft  but  with  morphing  airfoil  capability,  and  the  “geometry”  series 
represents  a  study  where  wing  area,  wingspan,  taper  ratio,  and  sweep  angle  are  free  to  vary  independently 
of  each  other.  Clearly,  a  wing  capable  of  telescoping,  chord  extension,  and  variable  sweep  behooves  a 
multi-role  platform. 
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Figure  4:  Comparison  of  Fixed-Wing  and  Morphing  Performances  [5] 

Even  a  brief  survey  of  aerodynamic  dividends  attainable  by  morphing  aircraft  promises  a  lucrative 
enterprise,  but  realizing  morphing  potential  is  a  multifarious  technical  challenge.  In  an  effort  to  resolve 
the  many  challenges,  the  first  two  phases  of  the  MAS  program  sought  to  achieve  four  technical  goals: 

1.  “Innovative,  active  wing  structures  that  change  shape; 

2.  Integration  and  aeronautical  use  of  advanced  sensors,  skin  and  structure  materials,  internal 
mechanisms  and  distributed  power  sources; 

3.  Advanced  capabilities  for  the  military  community; 

4.  Advanced  shape-changing  materials,  efficient  actuators  and  sophisticated,  smart  mechanisms.” 

[6] 

According  to  the  MAS  program  manager  Terrence  Weisshaar,  morphing  aircraft  wings  would  ideally 
have  at  least  a  50  percent  area  change,  thus  reconciling  the  reconnaissance  mission  requirement  of  a  wide 
wingspan  and  large  wing  area  with  the  combat  requirement  of  minimum  wing  area  for  speeds  of  Mach  2 
and  3.  As  the  MAS  program  entered  Phase  III,  two  contractor  teams,  Lockheed-Martin  Advanced 
Development  Programs  and  NextGen  Aeronautics,  sought  to  demonstrate  the  advantageousness  of 
morphing  by  comparing  performance  in  both  morphed  and  unmorphed  configurations  for  climbing  and 
turning  maneuvers.  Lockheed-Martin  developed  a  UAV,  dubbed  “Z-Wing”  (Figure  5),  with  two-position 
wings  that  essentially  fold  and  tuck  up  alongside  the  fuselage,  in  effect,  hiding  a  large  part  of  the 
planform.  Lockheed’s  1/8  inch  thick  skin,  fabricated  by  Shape  Memory  Polymer  (SMP)  technology,  was 
capable  of  stretching  100  percent  upon  stimulation  of  a  current  passing  through  it. 


3 

Approved  for  public  release;  distribution  unlimited. 


Figure  5:  Lockheed-Mar tin’s  Z-Wing  [6] 


Unlike  the  three-dimensional  morphing  of  the  Lockheed  UAV,  NextGen’s  morphing  concept,  called 
“Batwing”  (Figure  6),  consisted  of  in-plane  shape  change.  With  the  same  flexible  skin  challenge, 
NextGen  achieved  a  40  percent  area  change,  30  percent  span  change,  and  20°  change  in  sweep  while 
successfully  test  flying  a  remotely  piloted  vehicle  called  the  MFX-1.  In  light  of  the  technical 
competencies  hoped  to  be  gleaned  from  the  MAS  program,  progress  with  flexible  skin,  structural 
actuation,  and  suitable  flight  control  technologies  were  realized. 


Figure  6:  NextGen’s  Batwing  [6] 


1.2  Optimization 

As  stated  by  the  MAS  program,  the  technical  challenge  of  morphing  wing  design  will  involve  the 
multidisciplinary  integration  of  sensors,  actuators,  structures,  mechanisms,  and  skins;  a  challenge  far 
exceeding  simple  design  of  beams  or  four-bar  linkages.  Such  a  combinative  problem  is  likely  to  have 
intricate  solutions,  not  easily  approachable  on  the  basis  of  intuition  alone.  A  promising  design  approach  to 
a  multifaceted  problem  is  found  in  the  character  of  optimization.  Optimization  allows  a  designer  a  direct 
way  of  stating  any  number  of  objectives — quantities  to  be  minimized  or  maximized — and  design 
variables — parameters  that  can  assume  a  range  of  values  throughout  the  optimization  process.  Then,  once 
the  designer  has  identified  the  underlying  constraints  (usually  comprised  of  guiding  engineering 
principles  or  simple  geometric  relationships),  as  well  as  governing  bounds  on  the  design  variables  (set  by 
physical  or  other  reasonable  limitations),  optimization  solvers  iterate  through  a  set  of  equations  numerous 
times  until  a  stated  objective  converges  within  a  change  tolerance.  With  this  manner  of  problem 
formulation  and  assuming  the  mechanics  of  a  problem  have  been  correctly  formulated,  a  designer  focuses 
his  attention  on  identifying  the  most  desirable  objectives  to  meet.  Optimization  is  not  without  its 
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difficulties,  as  challenges  may  lay  in  interpreting  results  or  in  the  manufacturability  of  solutions,  but  it  is 
certainly  an  advantageous  tool  to  aid  the  morphing  challenge. 

Topology  optimization  seeks  the  optimal  layout  within  a  stated  design  domain  (Figure  7).  The  prescribed 
quantities  consist  only  of  the  loading  and  support  conditions,  the  fraction  of  the  design  space  to  be  filled 
with  material,  “volume  fraction,”  and  any  other  desired  design  restrictions  such  as  a  requirement  for  a 
hole  or  to  fill  the  boundary  of  the  domain.  The  shape  and  connectivity  are  not  prescribed  and  are 
determined  from  the  optimization.  The  topology  problem  is  called  a  distributed  parameter  system  because 
the  design  variables  represent  a  field  or  continuum  with  infinite  degrees  of  freedom.  Thus,  formulating  a 
calculable  problem  will  require  discretization  of  the  field  such  that  the  field  is  comprised  of  a  finite  set  of 
elements  each  with  a  finite  number  of  DOF.  The  selected  discretization  scheme  will  influence  the 
resulting  structure,  and  hence  incites  the  enterprise  of  obtaining  a  grid-independent  solution. 


Though  a  continuum  requires  discretization,  the  distributed  parameter  system  is  not  to  be  confused  with  a 
discrete  parameter  system  which  results  when  a  structural  problem  only  has  a  finite  number  of  design 
variables.  Such  problems  describe  naturally  discrete  systems  (i.e.,  the  sizing  of  truss  cross-sectional 
areas).  Topology  optimization  is  capable  of  driving  toward  a  discrete  or  a  continuous  solution  depending 
on  the  formulation,  and  either  solution  may  be  worthy  of  consideration,  depending  on  how  the  structure  is 
to  be  constructed  or  manufactured.  While  the  objective  function  may  encourage  a  discrete  or  a  continuous 
solution,  it  is  chiefly  seeking  to  minimize  a  particular  quantity.  What  that  quantity  represents  is  perhaps 
even  more  germane  to  the  structural  layout  than  is  the  discretization  scheme.  Of  the  many  objective 
functions  that  can  be  employed,  a  natural  starting  point  is  that  of  compliance  minimization,  which  is 
equivalent  to  stiffness  maximization.  Because  a  topological  approach  to  design  is  very  general  (i.e., 
determination  of  the  quantity  and/or  layout  of  beams,  trusses,  membranes,  actuators,  springs,  pivot 
locations,  etc.)  and  because  an  optimization  problem  can  be  formulated  by  an  assortment  of  objective 
statements  (i.e.,  compliance,  strain  or  potential  energy,  fundamental  eigenvalue,  buckling  load,  etc.), 
topology  optimization  is  highly  practicable  to  morphing  aircraft  structures. 

1.3  Problem  Description 

In  view  of  the  many  alluring  morphing  challenges,  the  structural  layout  of  morphing  MAV  wings 
subjected  to  a  gamut  of  aerodynamic  loads  as  the  MAV  executes  a  perching  maneuver  will  be 
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investigated.  A  perching  maneuver  is  a  common  landing  scheme  of  birds  in  which  they  dive  below  a 
landing  location,  exchange  kinetic  energy  for  potential  energy  while  pulling  up  to  the  perch,  and  then 
come  to  an  abrupt  vertical  landing  on  the  perch.  The  perching  maneuver  stimulates  much  interest  because 
it  is  a  miniature  of  a  multi-role  mission,  in  which  an  aircraft  leaves  a  reference  condition,  such  as  a  loiter 
wing  configuration,  “tucks”  its  wings  into  a  dive  configuration,  and  then  flares  them  back  out  both  to 
quickly  ascend  and  to  collect  momentum-reducing  drag. 

A  few  points  along  the  perching  trajectory  will  be  extracted  which  are  representative  of  the  extremes  of 
the  maneuver.  A  three-dimensional  estimation  of  the  aerodynamic  loads  will  be  made  and  applied  to  the 
morphing  wing  shapes.  By  constructing  a  compliance  optimization,  conceptual  wing  layouts  of  the 
structure  will  be  investigated  for  the  dissimilar  load  cases.  A  computational  tool  which  will  provide 
insight  into  the  design  mechanization  of  morphing  aircraft  structures  focusing  on  minimizing  the  energy 
required  to  achieve  the  desired  shape  change  was  also  developed.  Unlike  conventional  aircraft  design 
methodologies  which  aim  to  minimize  the  compliance  of  structures,  this  approach  will  leverage  the 
aeroelastic  effects  on  structures  to  adapt  the  configurations  in  an  effort  to  minimize  actuator  utilization. 
Recent  investigations  have  been  separately  conducted  for  mechanism  design,  topology  and  structural 
optimization,  actuator  placement,  and  energy  efficiency.  This  effort  aims  to  combine  these  disciplines  to 
produce  a  tool  for  efficient  integrations  of  actuators,  mechanisms,  and  elastic  structures  for  the  design  of 
morphing  aircraft. 
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Section  2.  METHODOLOGY 


This  chapter  details  the  theoretical  methodology  leading  to  the  development  of  the  MATLAB  codes  used 
to  predict  the  results  of  Section  1 .  The  initial  task  was  to  establish  a  geometry  generation  method,  and 
select  a  few  baseline  geometries.  Next  the  development  and  derivation  of  the  finite  element  models 
needed  to  solve  the  equilibrium  constraint  in  the  optimization  is  described.  A  perching  trajectory  with 
corresponding  flight  data  is  then  chosen,  and  the  process  for  estimating  the  resulting  aerodynamic  loads  to 
be  applied  to  the  FE  model  is  related.  The  compliance  minimization  objective  using  the  SIMP  model  is 
formulated,  and  the  OC  method  used  to  solve  the  optimization  problem  is  described.  Finally,  the 
optimization  methods  used  to  prescribe  the  optimum  placement  and  type  of  actuators  required  for  various 
shape  changes  is  presented. 

2 . 1  Geometry  Setup 

The  scope  of  the  current  investigation  is  limited  to  planform-changing,  two-dimensional  morphing,  akin 
to  NextGen’s  Batwing.  In  order  to  use  simple  bilinear  shape  functions  to  describe  the  geometry  of  the 
finite  elements,  planform  shapes  are  framed  by  straight  lines.  Additionally,  the  wing  is  comprised  of  any 
number  of  trapezoidal  sections.  This  restriction  ensures  that  the  leading  and  trailing  edges  have  the  same 
number  of  straight-edged  line  segments,  and  that  the  tip  chord  is  parallel  to  the  root  chord.  It  is  feasible 
for  any  geometry  that  meets  these  criteria  to  be  discretized  by  a  structured  mesh,  though  meeting  the 
criteria  does  not  guarantee  that  the  geometry  will  be  conducive  to  a  structured  mesh.  For  instance,  a  wing 
with  a  near-zero  taper  for  its  last  section  will  have  very  small  cells  at  the  tip,  which  may  lead  to  bad 
results  in  both  the  finite  element  analysis  and  the  aerodynamic  analysis. 

Figure  8  displays  a  MATLAB  generated;  three-sectioned,  general  wing  planform  mesh  in  terms  of  the 
script’s  input  variables.  The  code  provides  a  graphical  user  interface  (GUI)  for  easy  creation  of  a  wing  in 
any  number  of  configurations  (Figure  9).  More  details  of  the  generation  and  visualization  of  the  mesh  are 
presented  by  Elgersma  [7]. 
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Figure  9:  Wing  Geometry  GUI 

For  each  wing  configuration,  the  input  variables  define  a  root  chord  and  span  for  the  entire  wing,  and  a 
taper  ratio,  sweep  angle,  and  percent  of  the  total  span  for  each  wing  section.  The  mesh  is  created  by  linear 
interpolation  in  two  directions  with  the  number  of  grid  points  in  both  directions  specified  by  the  user. 
After  initial  development  of  the  script,  the  upper  profiles  of  a  NACA  4-series  airfoil  and  a  reflexed  airfoil 
were  hardcoded  into  the  script  to  add  camber  to  otherwise  flat  wings.  An  example  of  a  two-sectioned 
wing  with  camber  generated  by  the  GUI  and  is  shown  in  Figure  10. 
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Figure  10:  Example  Wing  Geometry 

In  order  to  capture  similar  planform  configurations  as  the  eagle  owl  (shown  in  Figure  1  through  Figure  3) 
is  capable  of  morphing  through,  a  rotating  mechanism  is  envisioned  that  can  sweep  about  5°  forward 
during  a  flare,  and  then  sweep  back  roughly  55°  to  60°  during  a  diving  dash.  The  trapezoidal  partitions 
and  generated  meshes  for  four  configurations  of  this  fictitious  geometry  (dubbed  “birdwing”)  are  shown 
in  Figure  1 1  through  Figure  14.  Table  1  gives  the  geometry  parameters  that  define  the  wing 
configurations.  In  later  analysis,  simple  rectangular  wings  with  and  without  sweep  are  also  used  for 
baseline  comparison  to  the  birdwing. 


-0.15 

o 


Figure  11:  Birdwing  in  Forward  Swept  Configuration  (m) 
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Figure  12:  Birdwing  in  Zero  Sweep  Configuration  (m) 


Figure  13:  Birdwing  in  Back  Swept  Configuration  (m) 
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Figure  14:  Birdwing  in  Dive  Configuration  (m) 
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Table  1:  Sectioned  Geometry  for  Multiple  Configuration  Birdwing 


Configuration 

5°  Forward 

0°  Sweep 

15°  Back 

Dive 

Chord 

m 

0.153 

0.153 

0.151 

0.136 

5 

Span (%) 

m 

0.069  (22.6) 

0.097  (31.8) 

0.070  (25.9) 

0.013  (9.7) 

a 

o 

Taper 

— 

0.988 

1.000 

1.040 

1.135 

V- > 

Sweep 

o 

-1.56 

0 

4.90 

54.74 

C3 

Oh 

QC  Sweep* 

o 

-1.17 

0 

3.68 

46.69 

Panels** 

# 

14 

19 

15 

6 

Chord 

m 

0.152 

0.153 

0.157 

0.154 

<N 

=tfc 

Span (%) 

m 

0.039(12.8) 

0.208  (68.2) 

0.200  (74.1) 

0.023  (16.9) 

a 

O 
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— 

1.014 

0.501 

0.483 

1.408 

V- > 

Sweep 

o 

4.97 

0 

-15.50 

54.74 

a 

Oh 

QC  Sweep 

o 

4.19 

5.24 

-9.95 

35.63 

Panels 

# 

8 

41 

45 

10 

Chord 

m 

0.154 

— 

— 

0.217 

m 

=tt 

Span (%) 

m 

0.197  (64.6) 

— 

— 

0.009  (6.7) 

a 

_o 

Taper 

— 

0.510 

— 

— 

1.081 

+- > 

t: 

Sweep 

o 

4.97 

— 

— 

30.89 

a 

Oh 

QC  Sweep 

o 

10.35 

— 

— 

5.97 

Panels 

# 

38 

— 

— 

4 
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m 

— 

— 

— 

0.234 

I 

Span  (%) 

m 

— 

— 

— 

0.011  (8.0) 

a 

Taper 

— 

— 

— 

— 

1.063 

V- > 
t; 

Sweep 

o 

— 

— 

— 

0 

a 

Oh 

QC  Sweep 

o 

— 

— 

— 

-18.99 

Panels 

# 

— 

— 

— 

5 

Chord 

m 

— 

— 

— 

0.249 

=tfc 

Span  (%) 

m 

— 

— 

— 

0.008  (6.2) 

eo 
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— 

— 

— 

— 

1.027 

V-> 

5_H 

Sweep 

o 

— 

— 

— 

-30.35 

a 

Oh 

QC  Sweep 

o 

— 

— 

— 

-38.07 

Panels 

# 

— 

— 

— 

4 

Chord 

m 

— 

— 

— 

0.256 

so 

=tt 

Span  (%) 

m 

— 

— 

— 

0.070  (52.5) 

a 

eO 

Taper 

— 

— 

— 

— 

0.040 

V-* 

5_H 

Sweep 

o 

— 

— 

— 

-72.65 

a 

Oh 

QC  Sweep 

o 

— 

— 

— 

-66.73 

Panels 

# 

— 

— 

— 

31 

Total  Span 

0.305 

0.305 

0.270 

0.133 

*QC  Sweep  is  the  quarter-chord  sweep 
**Panels  is  the  number  of  spanwise  panels 

2.2  Finite  Element  Model 

Previous  investigations  by  Joo  et  al.  [8]  [9]  and  Inoyama  et  al.  [10]  [11]  [12]  lacked  out-of-plane  loading 
capability  in  their  finite  element  models.  Therefore,  the  goal  of  this  section  is  to  develop  a  6  DOF  model 
that  can  support  membrane  (in-plane)  and  bending  (out-of-plane)  loads  and  be  easily  implemented  in 
MATLAB.  The  approach  undertaken  is  to  first  develop  membrane  (Section  2.2.2. 1)  and  bending  (Section 
2. 2.2. 2)  stiffness  matrices  independently,  and  then  superimpose  the  2  to  form  the  full  6  DOF  model 
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(Section  2. 2. 2. 3).  Section  2. 2. 2. 3  also  defines  a  transformation  matrix  so  that  non-planar  geometries  can 
be  analyzed.  After  the  element  stiffness  matrices  are  derived,  the  rest  of  the  finite  element  method  as  it 
relates  to  the  MATLAB  coding  process  is  described  in  Sections  2.2.3  through  2.2.5. 

In  the  following  sections,  the  finite  element  method  will  be  broken  down  into  the  following  eight  general 
steps: 

Step  1:  Discretize  the  Geometry 
Step  2:  Select  the  Element  Types 
Step  3:  Select  a  Displacement  Function 

Step  4:  Define  the  Strain/Displacement  and  Stress/Strain  (or  Equivalent)  Relations 
Step  5:  Derive  the  Element  Stiffness  Matrix 
Step  6:  Assemble  the  Global  Stiffness  Matrix 

Step  7:  Apply  Boundary  Conditions  and  Solve  for  the  Unknown  Degrees  of  Freedom 
Step  8:  Solve  for  the  Element  Stresses 

While  Steps  1  and  6-8  are  largely  the  same  for  the  different  element  types,  a  stiffness  matrix  must  be 
developed  for  each  element  type  in  Steps  2  through  5  (Section  2.2.2). 

2.2.1  Step  1:  Discretization 

Mentioned  in  the  discussion  of  the  geometry,  the  mesh  type  of  choice  is  a  structured  arrangement  of 
quadrilateral  elements  rather  than  an  unstructured  collocation  of  triangular  elements,  and  is  therefore 
easily  created  by  two-dimensional  linear  interpolation.  Another  advantage  of  a  structured  mesh  over  an 
unstructured  mesh  is  ease  of  coding.  Each  quadrilateral  element  has  four  definite  neighboring  cells  and 
can  be  mapped  to  a  rectangular  computational  domain,  allowing  a  programmer  to  simply  loop  through  the 
rows  and  columns  of  the  mesh.  The  elements  are  ordered  such  that  the  upper  left  corner  of  the 
computational  domain  is  the  first  element,  and  successive  elements  are  then  counted  down  through  the 
rows  and  then  right  through  the  columns  (see  Figure  18).  Using  general  quadrilaterals  rather  than 
rectangular  elements  allows  a  complex  geometry  to  be  fitted  with  a  structured  mesh.  However,  the 
disadvantage  is  that  cells  can  be  skewed  and  are  more  likely  to  have  large  aspect  ratios,  which  typically 
increases  the  inaccuracy  of  the  solution.  The  bird  wing  in  the  dive  configuration  (Figure  14)  demonstrates 
this  point  precisely,  where  highly  skewed  cells  dominate  the  sixth  wing  section  and  aspect  ratio  grows 
large  as  the  section  tapers  nearly  to  a  point. 

2.2.2  Element  Stiffness  Matrix 

2.2.2.1  Quadrilateral  Membrane  Element 

This  derivation  follows  the  nomenclature  and  methods  of  Reference  [13]. 

Step  2:  Select  the  Element  Type 

The  general  quadrilateral  membrane  element  considered  here  stretches  in  two  directions  and  thus  has  2 
DOF  for  each  of  its  4  nodes. 
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M= 


■;R}: 


(i) 


Writing  out  the  total  eight  displacements  gives 

[d}={ux  Vj  u2  v2  u3  v3  u4  v4}T  (2) 

Since  the  element  shape  is  an  arbitrary  quadrilateral,  a  transformation  from  the  physical  shape  to  a  square 
element  in  the  computational  domain  will  be  required  in  the  displacement  function  selection. 


Step  3:  Select  a  Displacement  Function 

In  an  isoparametric  membrane  formulation,  the  same  shape  functions  that  are  used  to  define  the  element 
shape  are  also  used  to  define  the  displacements  within  the  element.  Thus,  the  shape  functions  that  map  the 
natural  coordinates  r  and  s  of  any  point  in  a  square  element  to  the  global  coordinates  x  and  y  in  the 
general  quadrilateral  are  also  the  shape  functions  that  relate  the  internal  displacements  of  the  element  u 
and  v  to  its  nodal  displacements  {d}.  The  transformation  between  the  natural  coordinate  system  and  the 
global  is  depicted  in  Figure  15.  The  natural  coordinates  are  attached  to  the  element  and  rotate  with  the 
element.  They  need  not  be  parallel  with  the  global  coordinates  or  orthogonal  to  each  other.  The  four 
corners  and  edges  of  the  quadrilateral  are  bounded  by  +1  or  -1.  The  displacement  of  any  point  P  located 
at  (x,  y)  in  the  element  is  described  by  u  and  v. 


Figure  15:  General  Quadrilateral  Element  Mapping  to  Computational  Domain  [13] 


The  size  and  shape  of  the  quadrilateral  is  determined  by  the  eight  nodal  coordinates  x1?  y1?  x2,  y2,  x3,  y3, 
x4,  and  y4.  Therefore,  any  internal  coordinate  can  be  determined  by  bilinear  interpolation  of  the  natural 
coordinates  for  first-order  shape  functions.  The  global  coordinates  can  be  expressed  by  eight  unknown 
constants. 


x  =  Al  +  A2r  +  A3s  +  A4rs 
y  =  A5  +  A6r  +  A7s  +  Asrs 


(3) 
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The  A  coefficients  are  determined  by  evaluating  x  and  y  at  the  four  nodes  (-1,-1),  (1,-1),  (1,1),  and  (-1,1) 
in  natural  coordinates  and  (x^),  (x2,y2),  (x3,y3),  and  (x4,y4)  in  global  coordinates  and  simultaneously 
solving  for  the  A’s. 
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Solution  of  these  matrix  equations  for  the  A’s  and  substitution  of  the  expressions  for  the  A’s  back  into 
Eqs.  (3)  leads  to  the  following  mapping: 

x  =  [(l  —  r)(l  -  s)xx  +  (l  +  r)(l  -  s)x2  +  (l  +  r\  1  +  s)x3  +  (l  -  r)(l  +  s)x4  ] 

j;  =  ^[(1-rX1-‘sk  +  (i  +  rXi-A2  +(1  +  rX1  +  hy3  +(i-rX1  +  ^4] 

These  coordinates  can  be  expressed  in  matrix  form  in  the  following  manner: 
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1^4  J 

Here  the  shape  functions  \  relate  the  nodal  coordinates  ( x,.  y,)  to  any  (x,  y)  position  along  the  membrane. 
They  are  functions  of  the  natural  coordinates  (r,  s). 


v,  =f(l-rXl-s) 

V2=i(l  +  rXl-s) 

v3=f(i+<-Xi+h 

v<=i(i-<-Xi+h 

The  shape  function  partial  derivatives  with  respect  to  the  natural  coordinates  are 
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The  partial  derivatives  of  the  global  coordinates  then  become 
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(8) 
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\y*J 

The  displacement  functions  are  now  defined  in  the  same  manner  as  Eq.  (6)  with  the  shape  functions 
defined  by  Eqs.  (7): 


(9) 


(10) 


(11) 


where  u  and  v  are  displacements  parallel  to  the  global  x  and  y  coordinates.  Eq.  (11)  can  be  represented 
simply  as 


&'}  -  [N]{d}  (12) 

Step  4:  Define  the  Strain/Displacement  and  Stress/Strain  Relations 

Now  that  the  element  has  been  selected  and  the  displacement  functions  defined,  the  strain/displacement 
matrix  [B]  must  be  established  such  that 


(4  =  [B]{d] 


The  usual  relationship  between  strains  and  displacements  [13]  for  the  two-dimensional  case  is  given  as 


(13) 
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Eqs.  (14)  can  be  combined  into  the  following  matrix  equation: 
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(14) 


(15) 


or  more  succinctly 


(16) 

Since  the  global  coordinates  are  themselves  functions  of  the  natural  coordinates  (Eqs.  (6)  and  (7)),  the 
Jacobian  is  needed  to  express  the  partial  derivatives  of  the  global  coordinates. 
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Where 


(17) 


|  =  XrV*  ~  X«Vr 

Substituting  Eqs.  (17)  into  Eq.  (15)  yields  an  alternate  form  of  [D  ]: 
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Substitution  of  Eq.  (12)  into  Eq.  (16)  leads  to  the  definition  of  the  strain/displacement  matrix  [B]  by  Eq. 
(13): 

|B]  =  [01  [JV] 

(:)  x  8)  (3  x  2)  (2  x  8) 

Thus  [B]  is  computed  by  operation  of  [D  ]  on  [N]: 


[£(r,s)]  =  jh  [jjj  B,  B3  flj 


where  the  submatrices  of  [B]  are  given  by 


(-^f,r ) Vs  (  ^  t ,s  ) Ur  t) 

[S*\  =  0  U\»)*r  ~  (AYr)37* 

“  { ^i,r  )-T,5 

The  stresses  are  then  related  to  the  strains  by  the  constitutive  matrix  [D]  [13] 

n  -  im) 

where  [D]  is  defined  for  a  state  of  plane  stress  as 


[D}  = 


1  -if'1 


l  u  0 

v  1  0 

1  —  u 

(l  (i  - 


(20) 


(21) 


(22) 


(23) 


(24) 
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where  E  is  Young’s  modulus  and  v  is  Poisson’s  ratio.  Upon  substitution  of  Eq.  (13)  into  Eq.  (14),  the 
element  stresses  are  related  directly  to  the  nodal  displacements. 


W  -  [DPIM  (25) 

Step  5:  Derive  the  Element  Stiffness  Matrix 

To  derive  an  expression  for  the  stiffness  matrix  [k],  the  principle  of  minimum  potential  energy  is 
employed,  where  the  strain  energy  is  given  by 

u-rJJJ{4Ti^v  (26) 

V 

A  derivation  following  this  method  is  given  in  Reference  [13],  and  results  in  the  following  expression  for 
the  stiffness  matrix: 


M  =  JJJ\d\t\d\\dw 

V 

Since  [B]  does  not  vary  through  the  thickness  of  the  membrane  element,  the  stiffness  matrix  becomes  an 
area  integral. 


[k]  =  t  j j  \B\T\D][B)<l.yd>,  (2 

A 

However,  [B]  is  a  function  of  the  natural  coordinates  r  and  s,  and  not  the  global  coordinates  x  and  y  which 
define  the  volume  over  which  the  integration  is  to  be  performed.  Thus  the  variables  must  be  transformed 
using  the  determinant  of  the  Jacobian  matrix  J  (a  proof  for  a  general  transformation  is  given  in  [14]). 


fj  f(x,y)<hdy  =  JJ  J{r.s)\.l\drtlH  (29) 

A  A 

Since  the  area  is  bounded  by  ±1  in  the  natural  coordinates,  [k]  becomes 

r  =t  J1  J1  [B]T[D\[B]\J\dnh  (30) 

Since  [B]  is  a  rather  complicated  expression,  the  integration  determining  the  element  stiffness  matrix  is 
best  carried  out  numerically.  In  the  presented  work  the  Gaussian  quadrature  method  is  used  to  evaluate 
the  integral. 

2.2.2.2  Quadrilateral  Bending  Element 

A  general  quadrilateral  bending  element  formulation  was  not  found  among  the  finite  element  literature. 
Therefore  a  formulation  of  such  a  bending  element  was  derived  using  classic  Kirchoff  plate  theory  and 
using  the  same  shape  functions  that  defined  the  membrane  element  shape  in  the  previous  section.  The 
formulation  is  essentially  a  synthesis  of  Chapters  10  and  12  found  in  Reference  [13].  A  stand-alone 
derivation  of  the  bending  element  stiffness  matrix  is  provided  in  Reference  [7]  with  the  results  of  Steps  2 
through  5  repeated  here  for  completeness. 
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Step  2:  Select  the  Element  Type 

A  plate  element  can  translate  and  rotate  out-of-plane.  This  gives  the  plate  element  three  degrees  of 
freedom  per  node,  a  transverse  displacement  w  and  two  rotations  0X  and  0y. 


M  - 

(i, 

<h 

\  ;  {<*,}  =  < 

^ ^ 

dm 

k  ,l"  . 

A, 

Writing  out  the  12  displacements  gives 


(31) 


{</}  =  {  If,  0,,  tiy,  u-j  f)yj  irt„  Hr,,,  H,IW  «■„  Br„  J-  (32) 

As  with  the  membrane  element,  a  transformation  from  the  arbitrary  quadrilateral  to  a  square  element  in 
the  computational  domain  is  required. 


Step  3:  Select  a  Displacement  Function 

Although  the  same  shape  functions  that  are  used  to  define  the  membrane  element  shape  are  also  used  to 
define  the  bending  element  shape,  the  formulation  is  not  technically  isoparametric  because  the  nodal 
displacements  are  out-of-plane  and  have  their  own  shape  functions.  The  same  element  shape  mapping  is 
shown  in  Figure  16  for  the  bending  element,  where  representative  loads  are  shown  at  node  i. 


s 


i  >fd 


yt 


Figure  16:  Mapping  of  General  Quadrilateral  Element  to  Computational  Domain 


As  with  the  membrane  element,  the  size  and  shape  of  the  quadrilateral  is  determined  by  the  eight  nodal 
coordinates  by  Eqs.  (6)  and  (7),  and  obviously  have  the  same  partial  derivatives  (Eqs.  (9)  and  (10)).  Given 
that  there  are  12  degrees  of  freedom,  a  12-term  polynomial  is  selected  for  the  displacement  function  that 
enables  rigid  body  motion  and  constant  strain. 
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(33) 


The  rotations  follow  from  Eq.  (33)  since 


0X=^ 

% 


()lr 


e" _  a, 


r 


(34) 


The  displacements  can  be  put  into  the  form  of  Eq.  (12)  by  following  the  process  found  in  Reference  [7]. 


Step  4:  Define  Curvature/Displacement  and  Moment/Curvature  Relations 

Rather  than  relating  the  displacements  to  the  strains  (Eq.  (13))  which  vary  linearly  through  the  thickness 
of  the  bending  element,  displacements  are  instead  related  to  the  curvatures: 

{«}-[£]{<*}  (35) 

The  curvature/displacement  matrix  [B]  is  calculated  by  the  process  explained  in  Reference  [7].  Because 
the  stresses  also  vary  linearly  through  the  thickness  of  the  bending  element,  the  curvatures  are  related  to 
the  moments  instead. 


{M}  =  [niw 

where  the  constitutive  matrix  [D]  is  defined  as 


Here  D  is  expressed  as 


i  v  o 


[D\  =  D 


v  1 


0  0 


0 

1  -  V 
9 


D  = 


Et 3 

12(1  -  v2) 


Combining  Eqs.  (35)  and  (36)  relates  the  moments  to  the  displacements. 


(36) 


(37) 


(38) 


{M}  =  [D][0]{ri}  (39) 

Step  5:  Derive  the  Element  Stiffness  Matrix 

Step  5  must  be  modified  slightly  from  the  process  followed  in  the  membrane  derivation,  since  moments 
are  computed  rather  than  stresses.  The  potential  energy  of  the  plate  is  given  as 


u  = 


(  3  lj  A’j  f  \  1  ijK y  - \~  h  .m.,1  ) I ^ -4 


(40) 


Minimizing  the  potential  energy  leads  to  a  slight  difference  from  Eq.  (28),  where  thickness  no  longer  sits 
in  front  of  the  integral. 
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(41) 


[fc]  =  j j  \B]T[D}[B]dxdy 

A 

The  rest  of  Step  5  is  the  same  for  the  bending  element  and  results  in  a  [12  x  12]  stiffness  matrix. 

2.2.23  Combined  Quadrilateral  Membrane-Bending  Element 

A  shell  element  is  a  thin  curved  surface  with  six  degrees  of  freedom  at  each  node.  When  the  radius  of 
curvature  approaches  infinity,  the  shell  element  approaches  a  flat  plate.  When  a  surface  is  modeled  by 
many  small  shell  elements,  each  element  may  be  considered  as  a  flat  plate  with  a  particular  orientation  in 
space.  A  conventional  finite  plate  bending  element,  such  as  discussed  in  the  previous  section,  does  not 
allow  in-plane  membrane  deformation.  In  the  analysis  of  three-dimensional  structures,  membrane  and 
bending  stiffnesses  are  uncoupled  for  small  displacements.  Thus,  a  shell  element  can  be  formed  by 
superimposing  the  membrane  and  bending  stiffness  matrices  together  in  the  following  manner: 


[*]  = 
2U  *  20 


IM 

N  .  ^ 

i»i 

x  x  V> 

[0] 

*  12  %  8  | 

M 

12  x  12  - 

(42) 


where  subscripts  m  and  b  denote  membrane  and  bending  deformations  of  the  shell  element.  In  this 
formulation,  the  bending  and  membrane  elements  do  not  share  common  degrees  of  freedom  and  hence  are 
uncoupled.  Note  that  this  results  in  a  stiffness  matrix  of  order  [20  x  20],  or  five  degrees  of  freedom  per 
node.  The  element  structure  stiffness  equations  are 
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(43) 


For  a  single  shell  element,  there  is  no  twisting,  in-plane  rotation  0Z  (referred  to  as  a  “drilling”  rotation), 
but  for  an  assemblage  of  elements,  a  bending  rotation  in  one  element  may  result  in  a  drilling  rotation  in  an 
adjacent  element.  To  account  for  this,  the  drilling  degree  of  freedom  must  be  added  to  each  node,  raising 
the  order  of  the  stiffness  matrix  to  a  [24  x  24]. 
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(44) 


Since  there  is  no  stiffness  associated  with  the  drilling  degrees  of  freedom,  the  stiffness  matrix  will  be 
singular  if  all  the  elements  joining  a  node  are  coplanar.  If  flat  geometry  is  to  be  analyzed  by 
superposition,  the  on-diagonal  null  matrix  is  replaced  with  an  artificial  stiffness  matrix  which  causes  the 
twisting  rotations  to  produce  corresponding  moments  Mz .  One  suggested  fictitious  stiffness  matrix  [15]  is 
the  following: 
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(45) 
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where  a  is  a  small  number  such  as  0.3,  E  is  Young’s  modulus,  and  V  is  the  elemental  volume.  The 
element  volume  can  be  calculated  directly  from  its  coordinates  by  the  following  formula: 
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where  t  is  the  element  thickness.  The  current  order  of  the  displacement  matrix  i 
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Following  the  degrees-of-freedom  indexing  scheme  that  will  be  described  in  Section  2.2.3,  the 
displacement  vector  should  be  ordered  as  follows: 


(46) 
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Before  rearranging  the  combined  stiffness  matrix  [k],  the  membrane  and  bending  stiffness  matrices  can  be 
further  broken  down  into  submatrices. 
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The  fictitious  stiffness  matrix  of  Eq.  (45)  can  be  rewritten  slightly  as  follows: 
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(51) 
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where  kf  is  defined  as 


and  kf’  is  defined  as 
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Using  the  notation  of  Eq.  (49)  through  Eq.  (51),  the  total  membrane  and  bending  stiffness  matrix  is 
combined  [16]  in  the  following  manner: 
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(54) 


For  geometries  that  are  not  flat  with  entirely  coplanar  elements,  each  “shell”  element  has  a  unique 
orientation  with  a  local  axis  system,  indicated  by  a  *  over  the  variable.  A  transformation  matrix  [T]  can  be 
created  to  map  global  displacements  and  forces  to  corresponding  local  displacements  and  forces. 
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(55) 


{4  =  PM  ;  = 


The  local  and  global  stiffness  matrix  equations  are  given  as 


=  =  (56) 

Substituting  the  local  displacements  and  forces  of  Eqs.  (55)  into  the  local  stiffness  matrix  equation  yields 

=  (57) 

which  leads  to 

W-PT'WPM  (58) 

Thus  the  global  stiffness  matrix  is  computed  via  the  following  relation: 

1*1 -m  i\k]\T]  (59) 

Since  [T]  is  an  orthogonal  matrix,  the  transpose  of  [T]  can  be  used  instead  of  its  inverse  to  save 
computational  time. 

m  =  pf  torn  (60) 

Local  forces  can  also  be  mapped  back  to  the  global  forces  using  the  transpose. 
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The  transformation  matrix  [17]  relating  global  degrees  of  freedom  to  local  is  expressed  as 
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where  the  ly  ’s  are  the  direction  cosines  given  as 


(62) 
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Eq.  (62)  transforms  the  displacements  at  a  particular  node,  for  example,  at  node  i: 


«}  =  [T„K<M 

Thus,  the  transformation  for  a  four  node  element  becomes 


(64) 


24 

Approved  for  public  release;  distribution  unlimited. 


II 


[T]  = 


irn]  o  o 

G  x  G 

o  [rnj  o  o 

6  x  G 

o  o  [r„]  o 

G  x  G 

Q  <»  0  fr„] 

fi  X  l> 


The  direction  cosines  can  of  course  be  calculated  from  the  familiar  calculus  formula 
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However,  determining  the  unit  vectors  is  not  as  straight  forward.  The  global  unit  vectors  are  simply 


(66) 
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but  the  unit  vectors  of  the  natural  coordinates  will  require  some  effort  to  calculate.  The  element  normal 
unit  vector  et  can  be  computed  from  the  cross  product  of  two  vectors  lying  in  the  element,  for  example, 
vectors  extending  from  node  1  to  nodes  2  and  4  of  the  element  as  shown  in  Figure  17. 
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In  order  to  determine  the  unit  vectors  er  and  es,  the  midpoint  of  the  element  and  the  midpoints  of  the  sides 
r  =  1  and  s  =  1  are  needed  (shown  as  P5,  P6,  and  P7  in  Figure  17).  The  x  and  y  coordinates  of  these  three 
points  are  calculated  by  evaluating  Eqs.  (6)  and  (7)  at  (r,  s)  =  (0,  0),  (1,  0),  and  (0,  1).  The  plane  of  the 
element  can  be  defined  by  its  normal  vector  and  one  of  its  four  sets  of  nodal  coordinates: 


(68) 


f'fxt-T  —  rt  )  +  f’tyiV  ~  Vl )  +  '  f;(  Z  —  *J.)  =  0 


(69) 
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With  the  x  and  y  coordinates  obtained  from  the  shape  function  mapping,  Eq.  (69)  can  be  evaluated  to 
determine  the  z  coordinates. 


■  i  —  . —  [<r  tJ  [r  —  j  i )  !  t  hf(jf  —  ify  j  (70) 

&tz 

Thus,  unit  vectors  er  and  es  are  the  normalized  vectors  calculated  with  the  coordinates  of  points  5,  6,  and 
7. 


e* 


<Ar  ~  P L,  Pflff  -  ft,,  /l:  ”  ftz) 

(Ptx  -  /V  /%  “  ft*,  Pt.  -  Psz) 


(71) 


With  the  unit  vectors  er,  es,  and  et,  Eq.  (66)  can  be  applied  nine  times  to  compute  the  direction  cosines  of 
Eqs.  (63). 


2.2.3  Global  Stiffness  Matrix 

Step  6  of  the  finite  element  method  involves  assembling  the  local  stiffness  matrices  together  to  form  the 
global  stiffness  matrix  for  the  structure.  But  before  the  global  stiffness  matrix  can  be  compiled,  the 
elemental  stiffness  matrices  must  each  be  evaluated  by  Eq.  (30),  Eq.  (41),  or  Eq.  (54).  A  procedure  is 
constructed  to  loop  through  each  element  in  the  computational  domain.  Mentioned  in  Section  2.2.1,  the 
elements  are  indexed  by  counting  down  through  the  rows  and  then  to  the  right  through  the  columns.  The 
degrees  of  freedom  are  also  numbered  in  this  fashion  along  the  nodes.  An  example  two  DOF  membrane 
mesh  is  shown  in  Figure  18.  The  degrees  of  freedom  corresponding  to  each  element  are  then  ordered  in  a 
vector.  The  arrows  in  Figure  18  indicate  that  the  degrees  of  freedom  are  counted  counterclockwise 
starting  with  the  bottom  left  hand  comer  of  the  element  (see  also  Figure  15).  Indices  representing  the 
column  and  row  number  of  the  element’s  position  within  the  mesh  facilitate  selecting  the  corresponding 
degrees  of  freedom  for  the  element  based  on  its  number.  Along  with  the  degrees  of  freedom,  the  nodal 
coordinates  corresponding  to  each  element  will  also  need  to  be  assembled  to  compute  the  strain- 
displacement  matrix,  and  ultimately  evaluate  the  stiffness  matrix.  The  nodal  coordinates  define  the  shape 
of  an  element,  which  alone  differentiates  one  element’s  stiffness  over  another’s. 
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Figure  18:  Computational  Domain  Indexing  Example 

Once  all  of  the  elemental  stiffness  matrices  [ke]  have  been  evaluated,  the  direct  stiffness  method  [13] 
forms  the  global  stiffness  matrix  [K]  by  directly  assembling  individual  element  stiffness  matrices.  In  other 
words,  elemental  stiffness  matrices  are  superimposed  to  form  the  total  structure  stiffness  matrix. 


[AT  =  £[*«] 

e=l 


(72) 


2.2.4  Displacement  Solution 

In  Step  7  of  the  finite  element  method,  the  boundary  conditions  consisting  of  the  external  forces  and  any 
constrained  degrees  of  freedom  are  applied  to  the  global  structure  stiffness  equation.  In  a  manner  similar 
to  the  global  stiffness  matrix,  the  global  nodal  loads  are  obtained  by  lumping  the  body  forces,  distributed 
loads,  and  concentrated  nodal  loads  at  the  proper  nodes  into  a  column  vector.  Section  2.3.2  describes  how 
the  aerodynamic  load  calculated  for  the  bird  wing  is  divided  among  the  nodes. 

■v 

{F}  -£(/,}  (73) 

e=l 

Resolving  body  and  surface  forces  into  concentrated  nodal  loads  usually  involves  a  single  or  double 
integral.  An  example  of  integrating  a  surface  load  can  be  found  in  Reference  [13].  The  fixed  and  free 
degrees  of  freedom  are  specified  in  two  different  vectors.  For  instance,  one  boundary  condition  choice  for 
the  bird  wing  is  to  clamp  the  root  chord.  In  this  case,  the  user  would  specify  the  degrees  of  freedom  along 
the  root  to  be  fixed  and  all  other  degrees  of  freedom  as  free.  A  more  detailed  description  of  the  process 
and  associated  MATLAB  program  is  discussed  in  Reference  [7].  An  example  output  is  shown  in  Figure 
19  for  a  two-sectioned,  3 -DOF  model. 
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Figure  19:  Computational  Domain  Indexing  Example 


With  the  fixed  degrees  of  freedom  specified,  the  corresponding  displacements  can  be  directly  set  to  zero 
in  the  total  structure  displacement  matrix  {d}.  Thus  the  global  structure  stiffness  equations  can  be  solved. 


|F}  =  [K]{di  (74) 

It  is  important  to  note  that  {F}  and  {d}  may  represent  vectors  of  moments  and  rotations  rather  than  forces 
and  translations.  With  the  free  DOF  specified,  Eq.  (74)  can  be  elegantly  solved  in  MATLAB  by  the 
following  command: 


d^frccdofS)  1)  =  fT€6dof8)\F{fT€€dof&^  1)  (75) 

2.2.5  Element  Stresses 

The  final  step  of  the  method  consists  of  visualizing  the  results  of  the  displacement  solution.  After 
calculating  the  displacements  {d},  the  elemental  strains  and  stresses  are  computed  from  Eqs.  (13)  and 
(25),  respectively.  For  the  two-dimensional  state  of  stress,  three  independent  stresses  exist  and  are 
represented  by  the  column  vector 


/  \ 


* 


M  =  < 

The  principal  stresses  are  determined  from  Mohr’s  Circle  [18]. 


kTto 


(76) 
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(77) 


(78) 


The  von  Mises  effective  stress  ovm  [19]  is  defined  as  the  uniaxial  tensile  stress  that  would  create  the  same 
distortion  energy  as  is  created  by  the  actual  combination  of  applied  stresses,  and  thus  treats  multiaxial 


stresses  as  if  they  were  due  to  pure  tensile  loading.  Von  Mises  stress  for  the  two-dimensional  case  can  be 
calculated  with  the  principal  stress. 


(79) 


When  plotting  a  structure’s  deformed  state,  the  element  patches  can  be  colored  by  the  von  Mises  stress  for 


effective  visualization. 

Validation  of  the  membrane  and  bending  models  against  classical  analysis  are  detailed  in  Reference  [7]. 
Comparisons  with  the  individual  membrane  and  bending  models  are  also  made  for  the  combined 
membrane -bending  model. 

2.3  Aerodynamics 

2.3.1  Perching  Maneuvers 

In  the  interest  of  capturing  wing  loads  on  a  MAV  during  a  similar  maneuver,  perching  maneuver  data  was 
selected  from  Reference  [20].  Four  points  along  the  maneuver  were  selected  from  which  the  aircraft  and 
aerodynamic  data  were  extracted.  The  four  selected  points  are  marked  on  Figure  20  through  Figure  23, 
which  give  the  data  curves  of  altitude  versus  range,  range  versus  time,  angle  of  attack  versus  range,  pitch 
versus  range,  and  velocity  versus  range.  Table  2  provides  the  data  estimated  from  the  figures  at  each 
point. 

The  four  points  are  intended  to  represent  the  more  extreme  aerodynamic  conditions  of  the  maneuver. 

Point  1  is  at  the  beginning  of  the  trajectory  when  the  MAV  is  likely  in  the  15°  back  sweep  position 
(Figure  13).  The  dive  can  be  said  to  have  begun  since  the  near  -100%  angle  of  attack  and  pitch 
derivatives  indicate  the  vehicle  is  pitching  down  fast  from  its  cruising  or  loitering  condition.  At  Point  2, 
the  MAV  is  midway  through  the  dive  and  experiences  the  peak  velocity  of  the  maneuver.  The  MAV  starts 
to  pitch  up  again,  though  the  pitch  angle  is  still  negative.  The  wings  are  in  the  dive  configuration  (Figure 
14),  much  like  those  of  the  eagle  owl  in  Figure  2.  At  Point  3,  the  MAV  is  at  the  bottom  of  the  dive  and 
begins  ascending.  The  pitch  rate  is  at  a  maximum  due  to  how  fast  the  MAV  is  still  moving.  The  wings  are 
transitioning  from  the  dive  to  forward  sweep  configuration,  and  so  the  configuration  with  no  sweep 
(Figure  12)  is  selected  at  this  point.  Finally  at  Point  4,  the  MAV  has  covered  98.5  percent  of  the  range  and 
is  at  a  low  velocity  as  it  is  pulls  up  into  a  vertical  stance.  Consequently  the  angle  of  attack  is  post-stall  and 
the  angle  of  attack  rate  is  at  a  very  high  value  of  2167s.  The  wings  are  fully  flared  with  forward  sweep 
(Figure  1 1)  so  that  they  can  collect  as  much  drag  as  possible  to  bring  the  vehicle  to  a  halt. 
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Figure  23:  Perching  Maneuver:  Velocity  vs.  Range 
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Table  2:  Data  for  Selected  Points  along  Perching  Maneuver 


Selected  Location: 


Range 

m 

Altitude 

m 

Time 

s 

Velocity 

m/s 

AOA 

o 

AOA  Rate 

°/s 

Pitch  Rate 

°/s 

Point  1 

Point  2 

-20 

-12 

5 

3 

0 

0.92 

10 

10.41 

6 

3.75 

-98.92 

0 

-98.82 

0 

Point  3 

Point  4 

-6 

-0.31 

1.78 

4.80 

1.46 

2.44 

10.11 

2.19 

10 

50 

14.6 

216.2 

68.3 

26.4 

2.3.2  Force  Estimation 

With  the  geometry  and  flight  conditions  selected,  a  method  for  estimating  the  wing  loading  was 
established.  A  vortex  lattice  program  developed  in  MATLAB  called  “Tornado”  [21]  was  the  product  of 
the  Master’s  thesis  referred  to  in  Reference  [22].  Figure  24  shows  the  main  menu  of  Tornado. 

The  first  menu  item  allows  the  user  to  define  a  new  geometry  or  load  one  of  the  many  existing  wings.  A 
new  geometry  is  built  by  defining  any  number  of  wing  partitions  in  the  same  manner  as  Figure  8  with  the 
parameters  given  in  Table  2.  Thus,  the  vortex  panel  meshes  and  the  finite  element  meshes  were  created  to 
be  identical  to  one  another.  The  geometries  of  the  fictitious  sweeping  mechanism  shown  in  Figure  1 1 
through  Figure  14  were  created  in  this  manner  with  “flat  plate”  airfoils.  Next,  the  flight  conditions  of  the 
four  points  selected  in  Table  2  were  created  and  saved  in  menu  item  2.  The  lattice  was  then  generated 
with  menu  item  5,  where  the  fixed  wake,  standard  vortex  lattice  method  was  selected.  Figure  25  shows  a 
sample  geometry  output,  and  Figure  26  shows  the  vortex  panels  generated  for  the  sample  output. 


TORNADO  Version  135  Release  version 
build  2010  03  20  14:07  UTC 
Main  Menu 


Input  cp  e  r  at i ons . 

Il|  .  Aircraft  geometry  setup 

[2]  .  Flight  condition  setup 

13]  .  Change  rudder  setting 

14]  .  Move  reference  point 

Lattice  operations. 

[5].  Generate  lattice. 

Computation  operations. 

[€] .  Processor  access 

Post  processing  and  interactive  operations. 

[7] .  Post  processing F  Result /Plot  functions 

[3]  .  Keyboard  access 

Auxi 1 i a ry  ope  r at i ons . 

110]  .  About  /'  Release  Info 
[100] .  Help  files 
10] .  Exit  Tornado 

Please  enter  choice  from  above : 


Figure  24:  Main  Menu  of  Tornado  Software 
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Figure  25:  Example  of  Tornado  Geometry  Output  [21] 


In  light  of  the  several  processing  options  given  in  menu  item  6,  only  static  computations  at  the  specified 
flight  states  were  performed.  Tornado  outputs  five  different  MATLAB  structures:  geo,  state,  lattice,  ref, 
and  results.  The  geo  and  state  structures  simply  return  the  geometry  and  flight  conditions  defined  by  the 
user.  The  lattice  structure  returns  the  mesh  data.  The  ref  structure  is  not  useful,  since  the  mean 
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aerodynamic  chord  and  center  of  gravity  location  are  only  needed  for  full  aircraft  analysis.  The  pertinent 
data  is  reported  in  the  results  structure,  which  gives  three  column  vectors  of  the  x,  y,  and  z  force 
components  acting  on  each  panel.  Since  vortex  lattice  methods  are  inherently  inviscid,  an  alternative 
method  for  viscous  drag  estimation  was  made.  An  auxiliary  processing  option  was  added  to  the  current 
version  of  Tornado  to  predict  zero-lift,  flat-plate  drag,  which  gives  an  estimate  of  the  skin  friction  drag 
CDo.  For  the  perching  maneuver,  the  form  drag  due  to  flow  separation  as  the  MAV  exceeds  the  stall  angle 
of  attack  has  an  even  bigger  impact  on  the  wing  loading,  since  at  90°,  the  loading  is  almost  entirely  due  to 
drag.  (Although  it  should  be  noted  that  drag  loading  is  still  a  bending  load,  as  is  the  lift  loading  at  low 
angles  of  attack).  Therefore  to  estimate  the  form  drag  as  it  varies  with  angle  of  attack,  the  following 
equation  was  used: 

D  -  i S  [(  4-  ( 1  —  C  '/in )  sin  a]  (80) 

where  qoo  is  the  dynamic  freestream  pressure,  S  is  the  planform  area,  and  a  is  the  angle  of  attack.  The 
effect  of  Eq.  (80)  is  shown  in  Figure  27.  In  this  fashion,  the  drag  is  exactly  equal  to  the  zero-lift  drag 
prediction  at  0°  angle  of  attack.  At  90°  angle  of  attack  or  when  CDo  is  equal  to  one,  the  drag  is  equal  to 
qooS,  as  CD  equal  to  one  implies. 

To  apply  the  aerodynamic  forces  calculated  from  Tornado  to  the  finite  element  model,  the  single  resultant 
force  per  each  panel  calculated  by  the  vortex  lattice  method  needed  to  be  divided  onto  the  nodes  of  the 
given  element.  Since  the  locations  of  the  resultant  force  acting  on  each  panel  were  not  known,  the  forces 
were  assumed  to  act  through  the  centroid  of  the  element. 


Following  this  assumption,  one  fourth  of  each  of  the  three  components  of  the  force  vector  was  applied  to 
the  four  nodes  of  the  element.  Thus  at  each  node,  one  quarter  of  the  force  acting  through  each  of  the 
adjacent  elements  to  that  node  were  added  together.  To  distribute  the  estimated  drag  force  over  the  nodes, 
the  total  estimated  drag  was  first  divided  among  the  chords,  weighting  the  division  by  the  local  chord 
length.  Then,  the  component  of  the  total  drag  for  a  given  chord  was  distributed  evenly  among  the  nodes 
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along  that  chord.  In  this  manner,  the  estimated  drag  force  at  each  node  was  then  added  to  the  lift  and 
inviscid  drag  loads  calculated  by  Tornado. 

2.4  Optimization 

2.4.1  Compliance  Objective  Using  the  SIMP  Model 

In  structural  topology,  the  maximum  stiffness  or  rigidity  of  a  structure  is  often  sought  by  distributing  a 
limited  amount  of  material  throughout  the  design  domain.  In  this  pursuit,  a  common  restatement  of  the 
maximum  stiffness  objective  is  minimization  of  flexibility.  The  basic  relationship  of  flexibility  to  rigidity 
can  be  seen  with  the  fundamental  linear  spring  in  elementary  mechanics,  where  the  flexibility  is  simply 
the  reciprocal  of  the  spring  stiffness. 

In  order  to  minimize  the  flexibility  of  a  structure,  a  suitable  measure  of  flexibility  must  be  selected  to 
formulate  the  objective  of  an  optimization  problem.  Though  there  are  a  number  of  measures  of  flexibility, 
the  compliance  of  a  structure  will  be  used  here.  Compliance  can  be  thought  of  as  the  strain  energy  [18]  in 
a  deformed  solid  body.  Provided  no  energy  is  gained  or  lost  in  the  form  of  heat,  the  work  required  to 
deform  a  body  is  equivalent  to  the  strain  energy  of  that  body  in  its  deformed  state.  Work  is,  of  course,  the 
product  of  force  and  displacement.  Thus,  for  a  discrete  finite  element  formulation,  compliance  is  simply 
{F  }T  {d},  where  {d}  are  the  nodal  displacements  and  {F  }  are  the  external  forces  at  the  nodes  (see 
References  [23]  and  [24]  for  a  continuous  formulation  of  compliance).  Using  compliance  as  a  measure  of 
flexibility  has  the  advantage  that  it  is  a  convex  function  of  the  design  variables,  whereas  a  measure  such 
as  {d}T  {d}  is  a  nonconvex  function  of  the  design  variables.  Also,  a  body  where  compliance  has  been 
minimized  will  have  the  same  specific  strain  energy  in  each  element. 

Though  the  flexible  skin  process  summarized  in  Section  2.2  suggested  that  skin  design  should  include 
Young’s  modulus,  shear  modulus,  thickness,  and  possibly  Poisson’s  ratio  in  the  design  variables,  only 
element  thickness  is  considered  here  in  the  design  variables  as  a  simpler  starting  point.  Thus,  for  a  finite 
element  discretization  with  N  elements,  the  design  variables  are  given  in  the  vector 


P2  ■ - ■  Pe 


(81) 


where  pe  is  the  thickness  of  element  e.  It  will  be  discussed  later  whether  only  including  element  thickness 
in  the  design  variables  leads  to  results  that  can  be  correlated  to  skin  design. 


With  the  stated  objective  and  design  variables,  the  design  constraints  must  now  be  specified.  The  nodal 
displacements  {d}  in  the  objective  function  will  be  solved  by  the  finite  element  method  (Eq.  (75)),  and 
therefore  a  simultaneous  optimization  formulation  includes  Eq.  (74),  repeated  here,  as  a  constraint. 


in  =  AM 

Substituting  in  Eq.  (72)  in  for  the  global  stiffness  [K]  gives 


m 


i>,ij  w 


(82) 


(83) 
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As  stated  in  the  opening  sentence  of  this  section,  compliance  minimization  is  only  beneficial  if  a  limit  on 
the  material  distributed  throughout  the  design  domain  Cl  is  imposed.  In  the  absence  of  such  a  constraint, 
compliance  minimization  would  simply  lead  to  infinite  thickness  throughout  the  entire  design  domain;  but 
with  a  limited  volume  V  of  material,  the  optimization  must  “economize”  the  material  distribution. 
Expressing  the  material  constraint  as  an  inequality  gives: 


N 

V  P,ar  <  V  (84) 

e=l 

where  ae  is  the  area  of  element  e.  In  matrix  form,  this  constraint  is 

{p}T{«-}  -  V  <  0  (85) 

Additionally,  for  a  physical  solution  to  occur,  the  design  variables  must  be  confined  by  an  upper  and 
lower  bound. 


0  <  Pin  in  <  f\  <  Pinas*  1  =  1 . V  (86) 

Though  pmin  and  pmax  can  be  considered  inequality  constraints,  most  optimization  solution  methods  take 
into  account  optimization  variable  bounds  without  explicitly  adding  them  as  constraints.  Design  variable 
bounds  are  sometimes  referred  to  as  “box”  constraints. 


Putting  the  objective  function  and  constraints  together,  a  simultaneous  compliance  minimization  with  two 
constraints  is  given  by 


min 

p  d{p\ 


(HTW 


S.t. 


{<*}  =  m 


{py  {«}  -  v'  <  o.  n  <  pmi„  <  pe  <  pmaj. 


e=  1 . N 


(87) 


A  better  programmatic  formulation  of  the  problem  is  to  nest  the  equilibrium  constraint  in  the  compliance 
objective,  thus  making  the  objective  statement  a  function  in  the  design  variables  only,  and  not  in  d(p )  as 
well.  This  also  reduces  the  number  of  constraints  down  to  the  single  material  constraint. 


min  dip) 

p 


s.t.  {p}r{«}  -  V  <  0, 

where  the  compliance,  c,  is  given  as 


0  <  Pmin  <  Pt  < 


r|p)  =  { F} 1  {(/}.  wlitTc'  {</}  solves: 


L. 


m  =  {n 


(88) 


(89) 
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While  the  values  of  the  design  variables  are  allowed  to  range  between  pmin  and  pmax,  the  optimization 
problem  of  Eq.  (88)  is  really  a  sizing  problem,  as  defined  in  Section  1.2.  In  a  topology  design,  the  optimal 
placement  of  a  material  within  the  design  domain  Q  is  desired  such  that  points  in  space  will  either  have 
material  or  they  will  be  void.  In  this  case,  the  structural  geometry  is  rendered  as  a  black  and  white  image, 
albeit  pixelated  due  to  the  finite  element  discretization.  Thus  topology  optimization  is  seeking  to 
determine  the  optimal  subset  Qmat  of  material  points.  Ideally,  the  element  stiffness  matrices  would  be 
defined  piecewise  over  the  domain  Q  in  the  following  manner: 

{A  :j€  QnmJ 

(90) 

n  r  £ 

where  [kg(x)]  indicates  that  the  element  stiffness  matrix  varies  over  the  domain.  In  this  distributed, 
discrete  0-1  problem,  the  material  constraint  is  represented  as 


/  i,, . in  \\in ir . )  i  (9i) 

Jq 

Unfortunately,  a  number  of  problems  arise  with  the  character  of  a  discrete  problem;  allowing  finite 
elements  to  have  zero  stiffness  leads  to  singularities  in  the  solution  process.  Also,  optimization  algorithms 
where  design  variables  take  only  discrete  values  are  not  efficient  for  large-scale  problems  with  lots  of 
variables.  Thus,  a  model  that  allows  intermediate  values  in  the  design  variables  but  also  penalizes  them, 
the  so-called  Solid  Isotropic  Material  with  Penalization  (SIMP)  model,  is  used. 

In  the  SIMP  model,  the  constitutive  matrix  is  directly  penalized.  For  the  membrane  element,  the 
constitutive  matrix  (Eq.  (24))  takes  the  form 
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(92) 


where  the  penalty  q  is  some  constant  (typically  >1).  Since  the  upper  bound  pmax  is  set  to  1,  the  physical 
thickness  t  is  still  included  in  the  constitutive  matrix  such  that  when  an  element  is  fully  present  (p  =  1), 
the  actual  physical  thickness  is  scaled  by  t.  The  quantity  pq  can  be  thought  of  as  the  “effective”  thickness. 
To  avoid  numerical  problems  induced  by  zero  stiffness,  the  lower  bound  on  the  thickness  is  pmin  =  8  ~  0  (s 
=  0.001  is  typical).  Figure  28  shows  the  effect  the  choice  of  exponent  q  has  on  penalization.  When  q  =  1, 
the  stiffness  of  an  element  is  proportional  to  the  amount  of  material  designated  to  that  element,  and  hence 
the  thicker  an  element  is,  the  stiffer  it  becomes.  This  is  not  the  case  for  q  >  1.  Putting  some  values  to 
Figure  28,  at  half  the  available  material  (p  =  0.5),  the  effective  thickness  is  only  about  0.2  for  q  =  3.  In  this 
case,  the  effective  thickness  does  not  reach  one  half  until  a  p  of  about  0.8.  What  choice  of  exponent  q 
performs  the  best  is  a  matter  of  numerical  experimentation. 
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Figure  28:  SIMP  Penalization  of  Element  Thickness 


However,  Bendsoe  and  Sigmund  point  out  [23]  that  in  order  for  a  material  to  be  realizable  from  a  2-D 
SIMP  model,  q  must  satisfy  the  following  relationship  with  Poisson’s  ratio: 


(93) 


Thus,  for  a  material  with  v  =  1/3,  the  smallest  q  is  3.  In  fact,  q  =  3  has  been  found  to  be  a  good  choice  for 
the  exponent,  and  will  be  used  here  unless  otherwise  stated.  For  a  plate  bending  element,  the  constitutive 
matrix  (Eq.  (37))  reveals  that  the  element  stiffness  already  has  cubic  dependence  on  the  thickness.  This 
implies  that  the  optimal  design  naturally  prefers  to  achieve  either  the  upper  or  lower  bound  on  thickness, 
and  an  additional  penalty  need  not  be  imposed  (consult  Reference  [23]  for  more  information). 

Rather  than  keep  the  design  variables  in  the  constitutive  matrix,  the  global  stiffness  matrix  can  be  more 
readily  seen  as  a  function  of  the  design  variables  when  they  are  placed  outside  of  the  element  stiffness 


matrix: 


(94) 


With  the  amendment  to  the  compliance  objective  with  the  SIMP  model,  the  optimization  statement  is 
modified  slightly. 


min  dp! 

p 


(95) 


s.t.  {p}r{fi}  -  V  <  0,  £  <  pc  <1.  t  =  1 . N 


where  the  compliance  is  given  as 


{</}  =  {F}  (96) 
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Once  {d}  has  been  calculated  from  the  finite  element  method,  a  substitution  can  be  made  for  the  global 
force  vector  {F  }  such  that 


T 

{,/}  (97) 

The  total  compliance  can  be  computed  by  summing  up  the  compliance  of  each  element,  much  in  the  same 
way  that  the  direct  stiffness  method  adds  element  stiffnesses  and  forces  together  to  obtain  the  global 
stiffness  and  force  vectors.  Thus  Eq.  (97)  becomes 

N  :V 
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e=l  r=l 

The  Karush-Kuhn  Tucker  (KKT)  optimality  conditions  [25],  which  give  the  first-order  necessary 
conditions,  are  now  derived  for  the  compliance  minimization  of  Eq.  (95).  The  Lagrangian  function  is 
given  as 


c(p)  = 
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where  X  is  the  Lagrange  multiplier  for  the  single  volume  constraint,  and  s  is  the  slack  variable  for  the 
inequality  constraint.  A  formal  derivation  of  the  compliance  gradient  with  respect  to  the  design  variables 
is  given  in  Reference  [24],  where  both  a  direct  analytical  method  and  the  adjoint  analytical  method  are 
employed  to  attain  the  following: 
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(100) 


The  first  term  in  Eq.  (100)  is  only  included  when  the  global  force  vector  is  a  function  of  the  design 
variables,  for  example,  when  the  weight  of  the  structure  is  included  in  the  loading.  With  the  modification 
of  [ke]  by  the  SIMP  model  and  excluding  thickness  dependency  of  the  loading,  the  compliance  sensitivity 
reduces  to 


=  -wr,H}r[MTK} 
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The  gradient  condition  for  the  design  variables  thus  becomes 
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and  the  volume  constraint  is  recovered  from  the  gradient  condition 
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A  feasibility  check  on  the  inequality  requires  that 


*2  >  o. 
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(101) 


(102) 


(103) 


(104) 


The  switching  conditions  are  derived  from  the  derivative  of  the  Lagrangian  function  with  respect  to  the 
slack  variable. 


As  -  0  (105) 

Finally,  the  Lagrange  multiplier  is  required  to  be  nonnegative  since  it  is  associated  with  an  inequality. 

A  >  0  (106) 

The  KKT  conditions  produce  enough  equations  for  the  unknowns  and,  in  theory,  the  equations  could  be 
solved  simultaneously.  However,  the  nested  compliance  optimization  problem  has  been  shown  to  be 
convex  [24]  (although  the  simultaneous  formulation  is  not)  and  can  therefore  take  advantage  of  convex 
algorithms. 

2.4.2  Optimality  Criteria  Method 

Though  there  are  many  ways  to  solve  the  compliance  optimization  problem,  a  general  procedure  for 
solving  the  nested  formulation  using  first  order  algorithms  is  described  in  the  following  steps: 

1 .  Start  with  an  initial  design  p°  and  set  the  iteration  counter  k  =  0.  A  typical  initial  design  is  p°  = 
Vol(Qmat)/V  ol(O). 

2.  Perform  a  finite  element  analysis  for  the  current  design,  solving  for  the  displacement  vector  d(pk) 
from  the  equilibrium  condition  K(pk)d(pk)  =  F  . 

3.  For  the  current  iteration,  calculate  the  compliance  c(pk),  the  constraint  function  (pk)T  a  -  V  ,  and 
the  corresponding  sensitivities. 

4.  Evaluate  an  appropriate  update  scheme  (i.e.,  by  solving  an  explicit,  convex  approximation)  to 
give  a  new  design  pk+1. 

5.  Update  the  counter  k  =  k  +  1  and  return  to  step  2  unless  a  stopping  criteria  has  been  met. 

Step  4  provides  leeway  in  selecting  an  update  scheme  for  the  design  variables  and  indeed,  there  are  many, 
such  as  SLP  [25],  SQP  [25],  CONLIN  [24],  and  MMA  [26].  However,  the  Optimality  Criteria  (OC) 
method  is  a  historically  older,  classical  approach  that  can  be  seen  as  a  special  case  of  the  explicit  convex 
approximation  method,  and  has  turned  out  to  be  very  efficient  at  solving  topology  optimization  problems. 

A  more  descriptive  presentation  of  the  OC  method  is  given  in  Reference  [24],  but  only  a  brief  description 
is  presented  here.  The  OC  method  takes  advantage  of  the  fact  that  the  material  volume  is  a  monotonically 
decreasing  function  of  the  Lagrange  multiplier  X.  Hence,  a  stationary  point  is  achieved  when  the  volume 
constraint  is  satisfied.  The  update  of  the  design  variables  is  calculated  by 
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such  that  the  volume  constraint 
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is  satisfied,  where  pke+1  is  now  a  function  of  A.  The  update  scheme  can  be  written  a  little  more  compactly 
by  substituting  the  compliance  sensitivity  (Eq.  (101))  into  Eq.  (107). 
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(109) 
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In  this  manner,  each  thickness  is  updated  independently.  The  min  and  max  functions  are  used  so  that  the 
box  constraints  are  not  violated.  The  exponent  r|  acts  like  a  damping  knob  on  the  update  process  and  is 
chosen  by  experiment  (a  value  of  0.5  is  typical).  Note  that  the  design  has  converged  when  the  quantity  in 
parenthesis  is  equal  to  unity,  or  if 


<u{!  'K'VVdf 


U, 


=  A  =  constant*  for  all  e  —  1 *  A 


(110) 


This  quantity  is  twice  the  specific  strain  energy  of  an  element.  The  implication  of  Eq.  (1 10)  is  that  the  OC 
method  is  closely  related  to  a  fully  stressed  design,  where  all  elements  have  the  same  strain  energy,  and 
therefore  the  same  stress.  However,  with  the  penalization  method,  the  specific  strain  energy  is  only 
constant  for  elements  with  intermediate  densities,  since  it  is  lower  in  regions  with  p  =  pmin  and  higher  in 
regions  with  p  =  pmax.  The  design  update  does  however  add  material  to  elements  with  a  specific  strain 
energy  higher  than  X,  and  it  removes  material  from  elements  with  a  specific  strain  energy  lower  than  X, 
assuming  the  bounds  on  p  are  not  violated. 


A  bi-sectioning  routine  is  used  to  implement  Eqs.  (108)  and  (109).  First,  an  upper  and  a  lower  bound  on 
the  Lagrange  multiplier  are  selected,  for  example,  XuPP  =  100,  000  and  X\ow  =  0.  Next,  the  middle  point 
between  those  bounds  is  set  as  Xmid.  The  update  is  computed  with  Eq.  (109)  using  Xmid.  The  volume 
constraint  is  then  checked  with  the  updated  variables.  If  the  constraint  is  violated  with  excessive  material, 
then  the  lower  bound  is  set  to  Xmid,  and  if  the  constraint  is  violated  with  an  inadequate  amount  of  material, 
then  the  upper  bound  is  set  to  Xmid.  This  process  continues  until  the  value  of  X  is  found  within  a  certain 
tolerance  which  yields  a  design  update  that  satisfies  the  volume  constraint. 

An  example  output  is  shown  in  Figure  29  where  the  rectangular  design  domain  topology,  representing  the 
right  wing  of  an  aircraft,  has  been  optimized  for  an  in-plane  load  applied  at  the  upper,  right-hand  corner 
of  the  domain.  The  boundary  conditions  are  fixed  such  that  the  root  chord  is  clamped. 


Figure  29:  Example  “SIMP”  Output 
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As  a  final  note,  several  complications  can  arise  with  the  formulation  given  in  this  section.  A  common  one 
that  severely  inhibits  a  good  solution  is  that  of  the  so-called  checkerboard  pattern.  A  checkerboard  pattern 
occurs  when  the  design  function  alternates  between  0  and  1  and  produces  a  numerically  artificial  stiffness. 
This  problem  is  described  in  detail  in  Reference  [23].  Sigmund  devised  a  mesh-independency  filter  in  his 
99-line  code  for  the  compliance  sensitivities  that  also  eliminates  the  checkerboard  pattern.  Only  one 
additional  parameter  used  for  specifying  the  radius  of  the  filter  is  added  to  the  method.  A  radius  of  less 
than  1  indicates  the  filter  is  turned  off.  The  problem  in  Figure  29  is  repeated  again  in  Figure  30, 
displaying  the  erroneous  results  that  occur  when  the  checkerboard  filter  is  turned  off.  A  detailed 
investigation  of  the  influence  of  each  of  the  parameters  on  solution  convergence  is  discussed  in  Reference 

[7]- 


Figure  30:  Example  “SIMPm”  Output  with  No  Checkerboard  Filter. 

2,5  Substructure  and  Actuation  Optimization 

Morphing  structures  are  comprised  of  active  and  inactive  components  that  work  together  to  generate  the 
desired  shape  change  in  the  structure.  Active  components,  such  as  linear  and  rotational  actuators,  convert 
energy  to  motion  to  perturb  the  structure.  Inactive  components,  such  as  variously  compliant  beams  and 
joints,  guide  the  structure  towards  the  final  shape  and  provide  the  necessary  structural  rigidity.  Since  the 
objective  of  this  energy  based  design  of  the  substructure  and  actuation  system  is  to  minimize  the  amount 
of  actuation  energy  required  to  achieve  the  desired  shape  change  within  a  flight  environment,  the 
influence  of  external  forces  must  be  included  in  the  analysis  model.  This  chapter  focuses  on  the 
development  of  the  analysis  model  used  in  the  optimization  routine  and  provides  insight  into  element 
formulation,  model  kinematics,  and  external  force  estimation. 

In  order  to  properly  represent  the  design  space  for  a  morphing  structures  actuation  system,  an  innovative 
approach  to  finite  element  models  must  be  conceived.  This  model  must  facilitate  the  designation  of 
actuators,  rigid  and  compliant  beams,  rigid  and  compliant  joints,  and  even  voided  elements  within  the 
design  space.  The  formulation  of  the  model  must  be  general  enough  such  that  elements  can  take  the  form 
of  any  combination  of  active  or  inactive  components. 

The  model  kinematics  must  represent  the  motion  of  the  morphing  structure  as  accurately  as  possible. 
Simplifying  loads  due  to  actuators  as  static  loads  is  not  acceptable  since  the  static  loads  will  not  take  into 
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account  the  rotation  of  the  actuators  as  the  structure  changes  shape;  however,  a  dynamic  finite  element 
method  is  not  conducive  for  the  analysis  model  since  it  will  be  solved  many  times  during  the  optimization 
routine.  Therefore,  a  modified  static  finite  element  method  must  be  used  to  provide  a  better 
approximation  of  the  morphing  motion. 

To  model  the  external  forces  on  the  structure,  an  accurate  but  quickly  solvable  aerodynamic  code  must  be 
used.  For  this  reason,  a  vortex  lattice  method  (VLM)  is  used  to  quickly  calculate  the  external  aerodynamic 
forces  acting  on  the  structure.  These  forces  are  then  splined  to  the  structure  and  additional  information  is 
used  to  provide  a  static  aeroelastic  analysis. 

2.5.1  Elements 

The  mesh  and  connectivity  information  is  used  to  define  the  substructure  elements  within  the  model. 

These  elements  can  be  classified  as  either  line  elements,  joint  elements,  or  attachment  elements  as  shown 
in  Figure  3 1 . 


2. 5.1.1  Line  Elements 

The  foundation  of  all  line  element  configurations  is  the  classic  Euler-Bernoulli  beam  element.  Line 
elements  can  take  on  various  configurations  depending  on  the  values  of  their  associated  design  variables. 
Examples  of  these  configurations,  presented  in  Figure  32,  include  frames,  trusses,  inactive  telescopes, 
active  telescopes,  actuators,  and  voids. 
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•  Frame  elements  are  classic  Euler-Bernoulli  beam  elements  that  exhibit  both  axial  and  bending 
stiffness. 

•  Truss  elements  are  classic  elements  that  exhibit  only  axial  stiffness. 

•  Inactive  telescopic  elements  are  elements  that  exhibit  a  very  low  axial  stiffness  allowing  them  to 
extend  and  contract  when  external  loads  are  applied;  however,  they  also  exhibit  a  meaningful 
bending  stiffness. 

•  Active  telescopes  are  similar  to  inactive  telescopes  except  they  also  provide  an  actuator  force. 

•  Actuator  elements  exhibit  neither  axial  nor  bending  stiffness  and  provide  only  an  actuator  force. 

•  Void  elements  exhibit  neither  axial  nor  bending  stiffness  and  they  do  not  provide  an  actuator 
force. 

All  line  elements  are  subject  to  a  maximum  allowable  stroke  constraint  to  prohibit  them  from  extending 
or  contracting  an  unrealistic  amount.  Additionally,  actuator  type  elements  are  subjected  to  a  predefined 
maximum  force.  Formulation  details  of  these  configurations  are  presented  in  the  following  section. 

2. 5.1.2  Joint  Elements 

The  foundation  for  all  joint  elements  is  the  rotational  spring  element.  Joint  elements  can  vary  their 
rotational  stiffness  through  the  values  of  their  associated  design  variables.  Depending  on  the  stiffness  of 
the  spring,  the  joint  element  can  be  classified  as  revolute  joints,  compliant  joints,  or  semi-rigid  joints. 
Examples  of  these  configurations  are  presented  in  Figure  33. 


Figure  33:  Joint  Element  Configurations 


•  Revolute  joints  are  joints  elements  with  very  low  rotational  stiffness.  These  joints  are  hinges  and 
are  not  able  to  support  moments. 

•  Compliant  joints  are  joint  elements  with  moderate  rotational  stiffness  and  are  able  to  support 
moments  at  the  cost  of  high  rotational  strain.  These  joints  are  best  described  as  spring  hinges. 

•  Semi-rigid  joints  are  joint  elements  with  very  high  rotational  stiffness.  These  joints  are  essentially 
fixed. 

2. 5.1.3  Attachment  Elements 

Attachment  elements  are  basically  line  elements  at  the  boundary  between  the  wing  and  the  fuselage. 

These  elements  are  oriented  such  that  the  axial  direction  is  parallel  to  the  fuselage.  Constraints  on  the 
elements  degrees  of  freedom  are  set  such  that  the  element  is  not  able  to  penetrate  the  fuselage.  Further 
discussion  on  the  formulation  of  attachment  elements  is  detailed  in  the  following  section. 

2.5.2  Element  Formulation 

As  previously  mentioned,  in  order  to  properly  represent  the  design  space  for  morphing  structures,  an 
innovative  approach  to  finite  element  models  must  be  conceived.  This  single  element  or  a  multitude  of 
elements  must  be  able  to  accommodate  any  form  of  line  or  joint  configuration  within  the  design  space. 
This  section  focuses  on  the  formulation  of  the  element  used  in  the  analysis  model  which  is  not  only  able 
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to  accommodate  all  combinations  of  line  and  joint  configurations,  but  also  decreases  the  solution  time  of 
the  model. 


2. 5.2.1  Line  Element 

The  foundation  for  the  line  element  is  the  Euler-Bernoulli  beam  element.  This  element,  combined  with 
design  variables,  pn  and  pi2,  multiplied  to  the  axial  and  non-axial  portion  of  the  element  stiffness  matrix, 
respectively,  is  able  to  accommodate  all  of  the  various  line  element  configurations  (see  Eq.  (1 1 1)).  The 
design  variables  are  constrained  to  values  between  zero  and  one,  which  effectively  changes  the  stiffness 
of  that  element  to  a  percentage  of  the  total  stiffness.  For  instance,  if  the  design  variable  pn  has  a  very 
small  value  and  pn  has  a  value  of  one,  then  the  element  has  low  axial  stiffness,  but  high  bending  stiffness 
and  could  be  an  inactive  or  active  telescope. 
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2. 5.2.2  Joint  Element 

The  foundation  for  the  joint  element  is  simple  rotational  spring  element  which  solely  has  rotational 
stiffness.  This  element  is  combined  with  a  design  variable  much  like  the  line  element,  which  is 
constrained  to  values  between  zero  and  one  (see  Eq.  (112)). 
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2. 5.2.3  Condensed  Element 

What  makes  this  element  formulation  unique  and  favorable  to  this  type  of  analysis  is  how  the  individual 
line  and  joint  elements  are  condensed  into  a  single  element  stiffness  matrix.  This  process  is  formally 
known  as  Guyan  Reduction  [27]. 

Figure  34  represents  how  typical  analyses  using  beam  and  rotational  spring  elements  would  be  modeled. 
This  system  requires  four  different  nodes  to  fully  define  the  design  space.  Notice  the  axial  DOF  is  absent 
from  this  derivation  since  it  is  not  coupled  to  and  therefore  has  no  effect  on  the  non-axial  portions  of  the 
element  stiffness  matrix. 
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Figure  34:  Typical  Joint-Beam- Joint  Model 


Setting 


h  = 


E-t 


and 


k2  = 


the  stiffness  matrix  for  each  component  can  be  formulated: 
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Where,  ai  and  a2  are  functions  of  the  design  variables  associated  with  the  rotational  spring.  If,  equals 
zero,  then  the  joint  is  a  hinge  joint;  however,  if  oii  approaches  infinity,  then  the  joint  is  fixed  joint.  Details 
of  these  functions  are  discussed  later  in  this  section. 

From  equilibrium  requirements,  Fyl  equals  Fyi  and  Fy2  equals  Fyj,  and  from  compatibility  requirements,  Vi 
equals  v,  and  v2  equals  Vj,  the  stiffness  relationships  for  the  entire  system  can  be  assembled: 
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To  remove  0zi  and  0zj,  use  the  condensation  Eq.(l  17). 
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The  resulting  element  stiffness  for  the  Joint-Beam-Joint  configuration  (sans  axial  stiffness)  presented  in 
Eq.  (121)  is: 
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To  demonstrate  how  this  element  is  able  to  adapt  into  all  the  desired  configurations,  some  special  cases 
are  demonstrated  below. 

Rigid  Joint-  Stiff  Beam-Rigid  Joint  (a±  =  oo ;  pl2  =  1;  a2  =  °°) 
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ai  ~  Pfik joint max  (124) 

Please  note  pj  stands  for  “joint  design  variable”  and  is  not  index  notation.  The  purpose  of  defining  in 
this  manner  is  to  control  the  rotational  stiffness  through  a  design  variable  that  is  constrained  to  values 
between  zero  and  one.  If  pj  equals  zero,  then  equals  zero;  and  if  pj  equals  one,  then  ai  equals  kjoint  max. 
The  magnitude  selected  for  kjoint  max  was  1,000.  This  value  produced  a  stiffness  of  99.7  percent  of  a  true 
fixed  end.  The  Pinned  Joint-Stiff  Beam-Rigid  Joint  element  stiffness  matrix  is  provided  in  Eq.  (125). 
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In  addition  to  the  previously  derived  in-plane  element,  an  out-of-plane  element  of  a  similar  derivation  and 
a  membrane  element  were  added.  The  out-of-plane  element  supports  the  aerodynamic  loads  generated  by 
the  vortex  lattice  method  and  the  membrane  element  serves  to  provide  internal  resistance  within  the 
structure.  The  membrane  element  is  a  basic  four  node,  bilinear  quadrilateral  (Q4),  plane  isoperimetric 
element  that  is  unable  to  support  bending  loads.  The  program  is  configured  to  accept  user  a  specified  pre¬ 
strain  value.  Figure  35  demonstrates  how  the  pre-strain  is  applied. 


Figure  35:  Membrane  Pre-Strain 


The  pre-strain  is  achieved  through  the  patch  and  shrinkfaces  functions  in  MATLAB.  In  Figure  35,  the 
black  region  represents  the  original  dimensions  and  the  red  region  represents  the  “stretched”  dimensions. 
The  internal  load  associated  with  stretching  the  membrane  is  then  applied  to  the  respective  nodes  in  the 
element. 

2.5.3  Vortex  Lattice  Method  by  Ring  Vortex  Elements/Aeroelastic  Effects 

To  capture  the  aerodynamic  effects  on  the  structure  for  a  given  flight  configuration,  a  vortex  lattice 
method  (VFM)  based  on  ring  vortex  elements  as  described  by  Katz  and  Plotkin  [28]  was  added.  Results 
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from  the  VLM  were  compared  with  other  commercially  available  codes  such  as  Tornado  and  matched 
favorably.  Examples  of  Cp  plots  for  a  flared  wing  configuration  and  swept  wing  configuration  generated 
by  this  VLM  are  presented  in  Figure  36.  The  Cp  distributions  across  the  wings  indicate  typical  results. 


Figure  36:  Cp  Distributions 

During  the  VLM  subroutine,  the  aerodynamic  influence  coefficient  (AIC)  matrix  used  for  static 
aeroelastic  analyses  was  calculated.  This  matrix  was  used  to  calculate  an  aerodynamic  stiffness  penalty 
that  was  subtracted  from  the  global  stiffness  matrix.  This  process  provides  a  simple  static  aeroelastic 
result.  Figure  37  provides  an  example  of  the  deflections  provided  by  the  swept  wing  configuration  static 
aeroelastic  analysis. 
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Figure  37:  Static  Aeroelastic  Analysis 


2.5.4  Optimization  Background 

The  topology  optimization  relied  on  two  primary  tools  to  achieve  results.  The  first  is  the  Method  of 
Moving  Asymptotes  (MMA),  which  served  as  the  backbone  to  the  optimization  process.  The  second  was 
the  SIMP  method,  which  drives  the  solution  to  an  “off : ’  or  “on”  state.  The  SIMP  method  is  important  for 
obtaining  discrete  solutions. 

Method  of  Moving  Asymptotes  (MMA)  [23]  [26] 

The  Method  of  Moving  Asymptotes  (MMA),  developed  by  Svanberg,  is  a  method  for  non-linear 
programming  specifically  useful  in  topology  optimization  applications.  Non-linear  programming  refers  to 
the  process  of  solving  optimization  problem  formulations  where  the  objective  function  or  some  of  the 
constraints  are  non-linear.  Consider  a  structural  optimization  problem  in  standard  form. 
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where  ^>LV  1  is  the  objective  function,  ^'  A  l  A  — 11  are  the  constraints,  and  x 
vector  of  design  variables. 


is  the 


A  general  approach  to  solving  this  problem  formulation  is  to  generate  and  solve  a  series  of 
convex  approximating  sub-problems  in  an  iterative  manner: 
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Define  a  starting  point  a-1"1  ,  and  initialize  the  iteration  index  k  —  1 1  . 
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Given  an  iteration  point  x 
i  =  J . ni 


calculate 


and  the  gradients 


Generate  a  sub-problem  P*'  by  replacing  the  original  objective  function  and  constraints 
by  a  convex  approximation,  based  on  the  calculations  of  the  previous  step. 

Solve  the  sub-problem  P*1 . 

Set  the  solution  of  this  sub-problem  as  the  next  iteration  point  x[  1 J . 

Increment  k  and  got  to  set  I. 


Repeat  this  process  until  a  convergence  criterion  is  fulfilled,  or  when  the  user  is  satisfied  with  the  current 


solution  A'  .  Essentially, 
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is  obtained  by  a  linearization  of  A  in  variables  of  the  type 


U; 
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or 


'  V  depending  on  the  signs  of  the  derivatives  of  at  A'1*1 .  and  ^Js  are  parameters  that  are 
usually  updated  between  iterations  and  are  generally  refferred  to  as  “moving  asymptotes.”  The  process  for 
updating  these  parameters  will  be  discussed  later  in  this  section. 


In  order  to  generate  the  convex  approximating  sub-problem,  a  Taylor  Series  Expansion  is  taken 
from  the  original  problem  formulation  at  design  point  k  : 
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where, 
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Ignoring  higher-order-terms  and  transforming  the  independent  design  variables  in  the  form: 
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the  Taylor  Series  approximation  takes  the  form: 


(129) 
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Applying  Eq.  (133)  and  Eq.  (131)  to  Eq.  (130),  the  Taylor  Series  approximation  takes  the  form: 
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The  previous  derivation  only  took  into  account  the  upper  asymptotic  variables.  Applying  the  same 
procedure  to  the  lower  asymptotic  variables  provides  a  similar  approximating  function. 
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These  two  approximations  can  be  combined  into  a  single  approximation  where  the  upper 

—  >U  — ^<tJ 

asymptotic  approximation  is  used  when  1 ' A;  and  the  lower  asymptotic  is  used  when  P .  This 
single  approximation  is  shown  below. 
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Furthermore,  it  can  be  verified  that  •’*  is  a  first  order  approximation  of  ■<  at  A'  meaning, 
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)  and  df  ^/dxJ  =  dfjdx,  at  x=x«f 
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To  verify  the  convexity  of  this  sub-problem  approximation,  the  second  derivative  of  -A  with 
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Thus,  since  *  fi  —  and  —  ,  A  is  a  convex  function.  Looking  more  closely  at  Eq.  (145) 

/  (ft  F  J  (Jtji  (Jtj 

one  can  see  that  the  closer  /  and  J  are  to  '  v  ,  the  larger  the  second  derivative  becomes,  the 

more  curvature  is  given  to  the  approximation,  and  the  more  conservative  the  approximation  befits  the 

L  ^  If  (ft  i-  (ft 

original  problem.  Conversely,  if  i  and  ;  are  selected  far  away  from  J  ,  then  the 

approximating  functions  become  more  linear.  In  the  extreme  case  of  Ar  —  —  00  and  ~  00  , 

then  the  approximating  function  becomes  identical  to  the  linear  Taylor  Series  approximation 
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which  is  commonly  used  in  the  “sequence  of  linear  programs”  method. 

With  the  convex  approximating  function  defined,  the  following  sub-problem  formulation  is  given 
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where,  the  parameters  J  and  *  J  are  “move  limits,”  which  are  used  to  avoid  the  possibility  of  an 
unexpected  division  by  zero.  To  avoid  a  division  by  zero,  the  move  limits  should  be  chosen  such  that, 
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(149) 
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Methods  for  updating  the  moving  asymptote  parameters,  J  and  j  ,  will  now  be 
discussed.  A  simple  method  for  updating  the  moving  asymptote  parameters  is  to  let 
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where,  is  a  fixed  real  number  typically  less  than  one.  In  this  example  the  moving  asymptotes  are  not 
dependent  on  k  ,  which  means  they  are  “fixed  asymptotes”  rather  than  “moving.”  To  leverage  the 
complete  power  of  this  algorithm,  one  must  allow  the  moving  asymptotes  to  move  between  iterations. 

The  rules  governing  how  the  asymptote  parameters  should  move  are: 

1 .  If  the  process  tends  to  oscillate,  then  it  should  be  stabilized  by  moving  the  asymptote  closer  to  the 
current  iteration  point. 

2.  If  the  process  is  monotone  and  slow,  then  the  approximation  needs  to  be  relaxed  by  moving  the 
asymptotes  away  from  the  current  iteration  point. 

An  implementation  of  this  methodology  is  as  follow: 
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The  methods  of  adjusting  the  moving  asymptote  parameters  presented  in  this  section  are  those 
recommended  by  Svanberg.  There  are  numerous  alternatives  to  adjusting  these  parameters  and  further 
work  is  suggested  to  determine  the  best  practices. 

2.5.5  Problem  Formulation 

Properly  defining  the  problem  formulation  is  imperative  for  achieving  useful  results  in  an  optimization 
routine.  One  must  first  begin  by  developing  a  descriptive  statement  for  the  problem  which  describes  the 
overall  objective  and  requirements  to  be  met.  For  the  purposes  of  this  study,  a  design  set  of  structural 
members  and  linear  actuators  are  desired  which  morph  a  wing  from  one  shape  to  another  using  the  least 
amount  of  energy.  Next,  analysis  tools  are  identified  which  are  able  to  mathematically  model  the  design 
domain.  Finite  element  analysis  and  VLM  routines  provide  the  structural  and  aerodynamic  models 
necessary  to  define  this  design  space.  The  design  variables  must  then  be  chosen  that  define  the  design 
space  and  allow  for  a  feasible  optimization  solution.  The  design  variables  in  this  study  alter  the  stiffness 
characteristics  of  the  elements  as  previously  described.  The  final  steps  are  to  define  the  objective  function 
to  be  optimized  and  constraints  that  adequately  bound  the  problem. 

Many  variations  of  objective  functions  and  constraints  have  been  explored  during  this  effort  and  will  be 
expanded  upon  in  the  subsequent  thesis.  Outlined  in  this  report  is  the  current,  most  promising  objective 
function  and  constraint  set,  which  permit  the  optimization  of  multiple  wing  configurations. 

Topology  optimization  of  morphing  aircraft  structures  is  a  rather  difficult  problem  to  define.  The  design 
space  is  rather  open  and  problem  is  defined  by  a  very  large  number  of  design  variables.  In  addition, 
multiple  objectives  such  as  shape  matching  and  actuator  minimization  tend  to  conflict  when  the  solver  is 
trying  to  reach  a  global  minimum  value.  One  may  use  various  weighting  or  normalization  procedures; 
however  experience  has  shown  these  methods  are  difficult  to  balance.  A  two-step  optimization  has 
demonstrated  the  most  promise  in  producing  reliable  results. 

The  first  step  in  this  optimization  process  is  to  generate  a  set  of  design  variables  which  morph  the  wing 
into  all  configurations  without  considering  the  minimization  of  the  actuators.  This  shape  matching 
objective  drives  the  nodes  on  the  outer  mold  line  of  the  wing  to  pre-determined  locations  while  limiting 
the  stroke  of  each  line  element  and  the  overall  volume  fraction  of  the  wing. 

The  second  step  sets  the  shape  matching  criteria  as  a  constraint  along  with  the  stroke  and  volume  fraction 
constraints  but  focuses  on  minimizing  the  actuator  usage,  as  well  as  a  function  that  drives  the  elements  to 
a  0-1  solution. 

The  objective  functions  for  the  two  step  process  are  defined  below: 
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Step  1 


minimize: 


fo  =  DuTSft  -  U(ph)2 

subject  to: 

fCq  =  Kil  —  F  =  0  Static  EqttHibrum 

fm  =  —  Efnax  <  0  Stroke  Unit 

f+v  ~  E  P/ +  E  Pi  “  Vmax  <  0  Volume  Fraction  Limit 

i€U  i€L2 

where: 


(154) 


pmm  <  pi  <  10  i  6  n,  } 2,  LI,  L2,  B 
-1.0  <  pj  <  1.0  ;  e  A 


Step  2 

minimize: 

/o  =  p2  +  0i-(di^fi.)3] 

subject  to: 

feq  =  KU  -  F  =  0  Static  Equilibrum 

fm  =  E*,  -  El,ax  <  0  Stroke  Unit 

f+v  =  E  Pi  +  E  Pi  ~  Vntax  <  0  Volume  Fraction  Limit 

i€L1  i<=L2 

f shape  =  E(  U1"^1  -  U(p)i)2  -  shapeError^ux  Shape  Mulching 

where: 

Pmiti  <  Pi  <  1.0  i  €  n,  n,  LI,  LI,  B 
-1.0  <  pj  <  1.0  j€  A 


(155) 


The  second  objective  in  Step  2  is  novel  for  this  type  of  optimization  process  because  it  essentially  applies 
the  SIMP  method  twice  and  penalizes  values  that  are  neither  zero  nor  one  since  intermediate  values  will 
not  generate  a  minimum  value  for  that  particular  objective.  The  closer  the  design  variable  is  to  a  value  of 
one,  the  smaller  the  objective  gets.  The  traditional  SIMP  method  works  great  for  topology  optimization 
problems  that  aim  to  minimize  compliance  given  a  limit  on  the  amount  of  volume  fraction  allowed.  The 
SIMP  method  works  in  a  minimum  compliance  problem  formulation  because  values  that  are  neither  zero 
nor  one  do  not  effectively  utilize  the  potential  gain  in  element  stiffness  and  therefore  “waste”  valuable 
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volume  fraction  allotment.  For  topology  optimization  of  morphing  aircraft,  the  objective  function  as 
defined  without  this  added  objective  can  be  minimized  without  producing  a  0-1  solution.  Just  as  the 
volume  fraction  limit  drives  the  minimum  compliance  problem  to  0-1  solution  using  the  SIMP  method,  so 
does  this  additional  objective  in  topology  optimization  for  morphing  structures.  In  addition,  this 
additional  objective  drives  the  solution  to  favor  complete  beams  with  rigid  joints  such  that  only  necessary 
telescoping  members  and  re  volute  joints  exist. 

2.5.6  Sensitivity  Analysis 

Generally,  sensitivities  with  respect  to  the  design  variables  are  required  by  optimization  algorithms.  When 
the  objective  function  and  constraints  are  functions  of  only  the  design  variables,  these  sensitivities  are 
fairly  straightforward  to  derive.  For  functions  that  depend  on  the  displacements  as  well  as  the  design 
variables,  the  sensitivities  can  be  derived  by  the  chain-rule.  These  expressions  will  then  contain  the 
derivatives  of  the  displacement,  which  can  be  obtained  by  taking  the  derivative  of  the  static  equilibrium 
equation: 


KU  =  F 


(156) 


The  following  section  describes  how  the  sensitivity  of  the  shape  matching  objective  function  is 
calculated. 

In  topology  optimization,  a  moderate  number  of  constraints  are  typically  used.  Therefore,  it  is  generally 
better  to  calculate  the  sensitivities  using  the  adjoint  method  where  the  derivatives  of  the  displacements  are 
not  calculated  explicitly.  The  displacement  at  the  ith  degree  of  freedom  can  be  expressed  in  terms  of  the 
global  displacement  by: 


uf  =  ij  u 


(157) 


Where  It  is  the  vector  with  1  in  the  ith  and  0  in  all  other  components.  Since  KU  —  F  =  0,  this  function 
can  be  added  to  Eq.  (157)  without  altering  its  value  and  can  be  written  as: 

Uj  =  if  U  +  Xf  (KU  -  F)  (158) 


Where  A*  is  any  arbitrary,  but  fixed  real  vector.  Therefore  the  sensitivity  of  the  displacement  component 
can  be  written  as: 


dllj 

dpe 


!r“+Ar^u  +  AfK|^_A  r|£ 


dpe 


dpe 


dpe 


Which  can  be  rearranged  as: 


dUj 

dpe 


t  t  3LI  x  SK  t  dF 


(159) 


(160) 
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Since  the  derivative  of  the  global  displacement  vector  cannot  be  obtained  directly  and  since  A*  is  arbitrary, 
one  can  set  the  first  term  to  vanish  by  applying  the  following  condition: 


(If +\JK)  =  0 


(161) 


Which  can  be  rearranged  as: 


KAi  =  — Jj- 


(162) 


Where  Aj  can  be  solved  for  and  applied  to  Eq.  (160)  in  order  to  obtain  a  manageable  displacement 
sensitivity: 


—  —  Aj—U  -  aJ— 

dpe  dpe  dpe 


(163) 


Using  Eq.  (163),  the  gradient  of  the  multi-objective  function  can  now  be  expressed  as: 


d/o 

dpe 


-2  E  {U\arset  -  Ui)%A  e  e  J 1,  J2,  LI,  L2,  A 

i€T  pe 


(164) 


By  applying  Eq.  (163)  to  Eq.(164),  one  can  obtain  the  full  objective  function  sensitivity: 


-2WX  E ( u]arget  ~  Ut)\J $-U  eej  1 


d/o 

dpc 


-2W^{U\arset  eej  2 

ieT 

-21V,  E  ( U‘arget  -  u,)\j $-U  ee  LI 


-up  Aj^-U 


e  e  L2 


21Vi(/p^_1  E  (Ufarget  —  UpXfFg  e  e  A 


(165) 
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Section  3.  RESULTS 


This  chapter  explores  wing  structure  layouts  on  a  conceptual  level  for  several  wing  shapes  under  various 
loading  conditions.  First,  the  aerodynamic  results  from  Tornado  of  the  birdwing  geometry  (Figure  1 1 
through  Figure  14)  at  the  four  points  in  the  perching  trajectory  are  presented.  Then,  the  compliance 
minimization  runs  for  each  birdwing  layout  are  compared.  Finally,  a  brief  grid  independence  check  is 
provided. 

3,1  Aerodynamic  Results 

The  birdwing  of  Figure  1 1  through  Figure  14  with  zero  camber  (i.e.,  flat  plate)  maneuvering  through  the 
perching  trajectory  of  Figure  20  is  considered  here.  The  wing  is  in  the  back  swept  position  at  Point  1,  the 
dive  position  at  Point  2,  the  zero  sweep  position  at  Point  3,  and  the  forward  swept  position  at  Point  4. 

The  Tornado  calculations  of  the  birdwing  at  the  four  points  along  the  perching  trajectory  are  summarized 
in  Table  3.  The  velocity  and  angle  of  attack  have  been  included  for  reference. 

Table  3:  Aerodynamic  Data  for  Birdwing  along  Perching  Trajectory 


Point  1 

Point  2 

Point  3 

Point  4 

Vel. 

m/s 

10 

10.41 

10.11 

2.19 

AOA 

o 

6 

3.75 

10 

50 

Drag 

N 

0.0176 

0.0033 

0.0796 

0.0456 

Side 

N 

0.0061 

0.0075 

-0.0007 

-0.0087 

Lift 

N 

0.459 

0.112 

1.241 

0.196 

Fx 

N 

-0.0304 

-0.0040 

-0.1371 

-0.1208 

Fy 

N 

0.00610 

0.00749 

-0.00066-0.00871 

Fz 

N 

0.458 

0.112 

1.236 

0.161 

cL 

— 

0.220 

0.076 

0.512 

1.701 

CD 

— 

0.0085 

0.0022 

0.0328 

0.3958 

Cy 

— 

0.0029 

0.0051 

-0.0003 

-0.0757 

Re 

— 

90054 

137712 

91412 

19987 

CDo 

— 

0.0101 

0.0082 

0.0101 

0.0113 

Swet 

m2 

0.0681 

0.0444 

0.0775 

0.0785 

Dvis 

N 

0.2368 

0.1077 

0.4410 

0.0885 

Normal  N 

0.0248 

0.0070 

0.0766 

0.0678 

Axial 

N 

0.2355 

0.1075 

0.4343 

0.0569 

The  lift  forces  at  first  glance  seem  to  be  inordinately  small,  as  they  are  on  the  order  of  1  N.  However, 
since  the  wing  is  a  flat  plate,  an  estimate  of  the  lift  can  be  made  assuming  a  lift  slope  of  2n.  For  the  back 
swept  configuration  at  Point  1,  the  lift  coefficient  is  then 

c,  -  2nn  -  27r(G7r/lSO)  =  1 1.0580  ( 

The  lift  per  unit  span  becomes 
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V  =  ciq^c  =  0.6580 

If  the  wing  were  rectangular  with  uniform  lift  per  unit  span  across  the  span,  the  total  lift  force  for  the 
wing  would  be 


-(1.224-4  Jfcff/m3)(10  inf*) 


(0.153  m)  =  6.1630  N 


(167) 


L  =  L1  *t>  =  (6,163(1  AT}(0.270  m)  =  1.66 1  N  (ie 

The  lift  is  in  fact  on  the  order  of  1  N  for  such  a  small  wing  at  a  very  low  speed  (about  23  mph).  Since  the 
wing  is  tapered  and  the  lift  distribution  is  certainly  not  uniform  across  the  span,  the  resulting  lift  values 
are  necessarily  lower  than  predicted  by  the  simple  preceding  equations.  As  expected,  the  highest  lift 
occurs  at  Point  3  as  the  birdwing  begins  its  ascent,  and  the  lift  at  Point  2  is  the  lowest,  about  a  tenth  of  the 
value  at  Point  3. 

The  drag  calculated  by  Tornado  is  solely  the  induced  drag,  and  therefore  is  related  to  the  square  of  the  lift 
coefficient.  The  higher  drag  values  occur  in  the  ascending  portion  of  the  trajectory,  when  the  lift 
coefficient  is  highest.  The  fact  that  the  drag  at  Point  4  is  lower  than  at  Point  3  is  due  to  the  very  low 
velocity  of  2.19  m/s.  In  retrospect,  selecting  a  point  along  the  trajectory  between  Points  3  and  4  where  the 
velocity  is  only  somewhat  reduced  and  the  angle  of  attack  is  relatively  high  would  have  produced  a  higher 
drag  load.  The  drag  at  Point  2  is  smaller  than  the  drag  at  Point  3  by  a  factor  of  24,  due  to  the  reduced  wing 
area  and  the  smaller  angle  of  attack. 

The  distributions  of  the  pressure  coefficient  differences  between  the  upper  and  lower  surfaces  of  the  wing 
are  shown  in  Figure  38.  The  pressure  coefficient  differences  are  much  greater  for  the  forward  swept 
configuration  at  the  end  of  the  maneuver,  where  the  angle  of  attack  is  very  high.  Contrariwise,  the 
pressure  coefficient  differences  are  very  low  for  the  birdwing  as  it  dives,  because  it  is  slicing  through  the 
air  at  a  low  angle  of  attack.  In  each  of  the  four  plots,  the  pressure  coefficient  difference  is  about  four  times 
that  of  the  general  distribution  (the  yellow  areas)  over  the  wing. 
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Figure  38:  Cp  for  Birdwing  Geometries  along  Perching  Trajectory 


The  distribution  of  the  force  components  and  magnitudes  is  shown  for  each  configuration  of  the  birdwing 
in  Figure  39  through  Figure  42.  The  forces  shown  include  the  addition  of  the  viscous  drag  estimates  to  the 
body  forces.  The  viscous  drag  is  split  into  a  normal  and  axial  component  due  to  the  angle  of  attack.  The 
positive  x-axis  extends  towards  the  tip  of  the  wing  in  the  spanwise  direction,  and  the  y-axis  extends 
towards  the  leading  edge  in  the  chordwise  direction. 
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(a)Fx  (b)Fy 


(C)  Fz  (d)  Fmag 


Figure  39:  Force  [N]  for  Birdwing  in  Swept  Back  Configuration  at  Point  1 

At  each  point  in  the  maneuver,  the  side  forces  (Fx)  are  relatively  negligible  and  not  expected  to  influence 
the  structural  layout  much.  The  axial  forces  stimulate  interest  because,  at  the  leading  edge,  the  forces  act 
in  the  positive  y-direction,  whereas  everywhere  else  on  the  wing,  the  axial  forces  act  in  the  negative  y- 
direction  (as  intuition  suggests).  The  fact  that  the  leading  edge  is  seemingly  being  pulled  forward  is  due  to 
the  domination  of  the  axial  component  of  the  lift  over  the  axial  component  of  the  drag.  The  lift  calculated 
by  Tornado  at  each  panel  over  the  surface  of  the  whole  wing  dominated  the  drag  in  the  axial  direction; 
however,  the  added  viscous  drag  reversed  the  net  force  over  the  all  of  the  wing  except  the  leading  edge. 
Hence  the  structure  is  being  stretched  in  both  axial  directions. 
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Figure  40:  Force  [N]  for  Birdwing  in  Dive  Configuration  at  Point  2 


The  loads  for  the  dive  configuration  are  the  on  the  order  of  10~4,  whereas  the  loads  for  the  other  cases  are 
on  the  order  of  1CT3;  thus,  the  dive  configuration  is  probably  not  a  good  design  point  since  ultimately  the 
loads  at  each  configuration  will  have  to  be  withstood  by  a  single  structure. 
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Figure  41:  Force  [N]  for  Birdwing  in  Zero  Sweep  Configuration  at  Point  3 

Clearly  the  load  on  each  of  the  configurations  is  leading  edge  dominated.  Thus  a  good  baseline 
comparison  for  resulting  structures  could  be  constructed  by  simply  distributing  a  load  over  the  leading 
edge  of  a  similar  geometry.  Also,  resultant  magnitude  contours  of  each  point  in  the  maneuver,  except  the 
last,  match  the  contours  of  the  normal  component  Fz ,  indicating  that  the  bending  loads  are  the  dominant 
loads.  The  flaring  of  the  wings  at  the  end  of  the  maneuver  should  be  bending  dominated  as  well,  but  Point 
4  that  was  selected  occurs  when  only  a  fraction  of  the  peak  velocity  remains.  Therefore,  the 
corresponding  drag  loads  are  low,  and  the  bending  dominance  is  not  captured.  In  all  cases  (except  perhaps 
in  the  case  of  the  dive),  the  trailing  edge  of  the  tip  of  the  wing  experiences  essentially  no  load  and 
suggests  that  corresponding  structures  may  clip  the  lower  tip  of  the  wing. 
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Figure  42:  Force  [N]  for  Birdwing  in  Forward  Swept  Configuration  at  Point  4 


Though  the  magnitudes  of  the  body  forces  of  the  birdwing  in  the  perching  trajectory  is  very  low,  the 
relative  magnitudes  are  of  greater  importance  to  the  optimization  process.  The  majority  of  the  relative 
magnitudes  of  lift  and  drag  seem  acceptable,  but  the  most  unreasonably  captured  aerodynamic  extreme  is 
the  drag  at  the  end  of  the  maneuver  which  should  be  a  bending-dominated  load.  The  beginning  of  the 
maneuver  is  also  dominated  by  bending  loads,  but  due  to  the  lift  and  not  the  drag.  If  the  optimization 
process  determines  any  hybrid  membrane -bending  structures,  it  is  anticipated  that  the  Point  3  flight 
condition  would  produce  this  effect. 


3.2  Optimized  Static  Structural  Layouts 

For  each  point  along  the  perching  maneuver,  three  sets  of  loading  scenarios  were  performed.  The  first  set 
consisted  of  only  the  membrane  forces  (Fx  and  Fy),  the  second  set  considered  the  bending  forces  (Fz ) 
only,  and  the  third  set  was  the  combined  membrane -bending  loading  scenarios.  Each  set  consisted  of  four 
cases  with  volume  fractions  ranging  from  0.2  to  0.5.  Unless  convergence  issues  arose,  each  case  had  a 
penalization  of  3,  a  filter  radius  of  1 .5,  and  was  run  with  a  move  limit  of  0.2  and  a  damping  of  0.5.  All  the 
runs  were  minimized  until  the  solution  converged  within  a  0.01  maximum  thickness  change. 

The  results  of  the  compliance  minimization  for  the  birdwing  at  Point  1  in  the  perching  maneuver  is  shown 
in  Figure  43  through  Figure  45.  Considering  first  the  purely  membrane  solution  of  Figure  43,  the 
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predominate  structural  features  are  the  two  members  extending  from  the  root  of  the  wing  in  the  form  of  a 
wishbone,  and  at  the  tip  of  the  wing,  a  second  wishbone  forms.  Even  though  the  predominate  membrane 
forces  occur  along  the  leading  edge,  the  leading  edge  is  not  built  up  across  the  whole  span  like  it  is  for  the 
topologies  of  Section  3.1.  This  is  due  to  the  force  along  the  remainder  of  the  wing  pulling  the  wing  in  the 
opposing  negative  chordwise  direction.  The  “wishbone”  is  a  natural  structure  to  resist  a  spreading  motion 
(the  converse  motion  of  squeezing  a  pair  of  pliers).  It  can  be  thought  of  as  antipodal  to  two  joined 
buttresses.  Had  all  of  the  force  in  the  chordwise  direction  been  directed  towards  the  trailing  edge,  the 
structures  would  be  expected  to  have  a  fully  supported  leading  edge.  The  second  main  members  are  those 
along  the  leading  and  trailing  edges.  As  evidence  that  the  two  wishbones  and  the  leading  and  trailing  edge 
members  comprise  the  main  structure,  each  time  material  is  added  to  the  wing  domain,  these  members  are 
thickened. 


(c)  Volf  =0.4  (d)  Volf  =0.5 

Figure  43:  Membrane  Structure  for  Birdwing  at  Point  1 

The  wishbone  structure  and  the  intermediate  members  form  a  scissor-like  structure,  and  consequently  the 
voided  areas  are  predominately  quadrilateral  in  shape.  This  is  in  contrast  to  the  usual  triangular  voids  in 
trusses. 

Another  striking  feature  of  the  structure  is  the  straight  protrusions  of  intermediate  thickness.  These 
battens  or  rods  point  towards  the  center  of  a  voided  area,  and  towards  the  unsupported  perimeter  of  the 
wing  domain.  Thus  they  are  comparable  to  spokes  on  a  wheel.  This  feature  is  a  result  of  a  distributed  load 
covering  the  en-tire  surface  of  the  wing,  and  is  absent  in  an  over-simplified  point  load  model.  At  a  fifty 
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percent  volume  fraction,  the  topology  finally  adds  new  members  in  the  form  of  another  smaller  wishbone 
inscribed  within  the  main  wishbone. 


(c)  Volf  =  0.4 


(d)  Volf  =  0.5 


Figure  44:  Bending  Structure  for  Birdwing  at  Point  1 


The  bending  structures  of  Figure  44  reveal  that  the  basic  form  of  the  bending-resistant  structure  is  not  a 
truss,  but  rather  a  conglomeration  in  the  form  of  a  beam.  The  importance  of  exploring  a  range  of  volume 
fractions  for  a  given  wing  shape  and  loading  scenario  is  evident  in  the  progression  from  Figure  44(a)  to 
Figure  44(d).  The  impression  of  the  50  percent  volume  fraction  case  is  that  the  bending  structure 
primarily  supports  the  leading  edge,  where  the  bending  forces  are  highest.  But  Figure  44(a)  shows  this  is 
clearly  not  the  case.  The  most  basic  or  crucial  structure  is  a  beam  that  is  situated  near  the  quarter  chord  of 
the  wing.  For  the  flat  plate  birdwing,  the  quarter  chord  is  also  approximately  the  center  of  pressure  and 
the  aerodynamic  center  of  the  wing,  where  the  aerodynamic  moment  is  zero.  As  the  volume  fraction 
increases,  material  is  added  along  the  sides  of  the  original  beam  of  Figure  44(a).  The  developing  shape  of 
the  beam  is  similar  to  the  shape  of  a  wishbone  (i.e.,  like  A-arms  in  a  suspension  system),  or  in  a  more 
crude  sense,  the  shape  of  a  delta. 


Like  the  rods  or  battens  in  the  membrane  structure,  the  bending  structure  also  develops  a  secondary 
feature  captured  by  the  intermediate  thicknesses.  Rather  than  straight  rods  with  a  regular  arrangement, 
branches  stem  from  the  central  beam  structure  in  arbitrary  directions,  and  then  further  divide  into  smaller 
networks  that  attempt  to  cover  as  extensive  a  region  as  possible.  The  branches  resemble  veins  in  the 
wings  of  insects,  and  as  the  high  density  of  nerves  in  bats  attest,  are  a  very  effective  method  of  stiffening  a 
membrane  skin. 
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The  combined  membrane  and  bending  structure  for  the  first  point  in  the  maneuver  is  shown  in  Figure  45, 
revealing  an  identical  structure  to  the  membrane  case.  However,  when  the  viscous  drag  is  removed,  the 
result  is  that  of  Figure  46,  which  reveals  a  hybrid  solution  of  the  two  isolated  membrane  and  bending 
cases. 


Figure  45:  Combined  Structure  for  Birdwing  at  Point  1 
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Figure  46:  Combined  Structure  for  Birdwing  at  Point  1  without  Viscous  Drag 

This  is  not  immediately  evident  in  Figure  46(a),  but  the  higher  volume  fraction  solutions  all  have  a  delta¬ 
shaped  “beam”  on  the  front  half  of  the  wing.  All  of  the  solutions  have  supporting  truss  structure  along  the 
trailing  edge  of  the  wing  that  then  connects  to  the  leading  edge  near  the  tip  of  the  wing.  As  material  is 
added,  the  general  topology  remains  unchanged,  however  the  beam  portion  of  the  structure  shifts  in  the 
chordwise  direction,  such  that  at  fifty  percent  volume  fraction,  it  is  centered  mid-chord.  Both  the  scissor 
pattern  from  the  membrane  solution  and  the  beam  extending  most  of  the  way  in  the  spanwise  direction  are 
attenuated,  and  the  emerging  structure  better  resembles  a  perimetric  truss  structure  with  spokes  protruding 
inward  to  a  hub  covering  the  upper  left  corner  of  the  wing.  The  vein-like  stiffeners  are  absent,  and  only  a 
few  rods  are  present.  While  the  hybrid  solution  is  interesting,  in  actuality  a  completely  membrane 
structure  supports  the  combined  loading  according  to  the  compliance  results.  The  larger  lifting  loads  and 
smaller  drag  loads  at  the  relatively  small  angle  of  attack  of  Point  1  would  suggest  a  more  bending- 
dominated  solution,  but  the  effect  of  the  distributed  viscous  drag  load  drives  the  solution  towards  a 
membrane  structure. 

As  previously  mentioned,  the  diving  portion  of  the  perching  maneuver  is  perhaps  not  a  good  design  point 
when  the  magnitude  of  the  loads  experienced  throughout  the  entire  maneuver  is  considered.  However,  the 
wing  shape  is  quite  distinct  from  the  other  three  configurations  and  thus  yields  unique  structures  worthy 
of  investigation.  The  resulting  structures  are  presented  in  Figure  47  through  Figure  49. 
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(a)  Volf  =  0.2  (b)  Volf  =  0.3 


(c)  Volf  =0.4  (d)  Volf  =0.5 


Figure  47:  Membrane  Structure  for  Birdwing  at  Point  2 

If  obscured  by  other  truss  members  in  the  back  swept  wing  configuration,  there  is  no  overlooking  the 
prominent  wishbone  of  Figure  47.  In  fact,  it  is  the  only  truss  structures  present  in  Figure  47(a)  and  (b).  In 
Figure  47(c),  a  smaller  wishbone  facing  opposite  of  the  main  wishbone  is  inscribed  within  the  main 
wishbone  again  producing  a  quadrilateral  void  and  a  scissor-like  structure.  Aside  from  this  small  amount 
of  “truss”  structure,  the  remainder  of  the  structure  consists  of  very  long  and  definitive  battens,  in 
comparison  with  those  in  the  back  swept  wing  structure.  In  fact,  the  battens  are  not  of  intermediate  value, 
but  are  maximum  thickness.  As  shown  in  Figure  40(b),  the  maximum  axial  load  occurs  at  the  leading 
edge  near  the  tip  of  the  wing,  and  accordingly  a  batten  extends  to  the  tip  of  the  wing.  The  region  is  not 
just  darkened  because  of  the  overly  refined  mesh.  In  light  of  the  poor  finite  element  mesh  of  the  dive  wing 
shape,  an  instructive  structure  develops. 

The  compliance  minimization  of  the  bending  structures  shown  in  Figure  48  was  troublesome,  and  only 
the  forty  and  fifty  percent  volume  fraction  runs  yielded  results  that  were  not  explicitly  erroneous.  The 
twenty  and  thirty  percent  volume  fraction  cases  were  never  close  to  converging,  even  for  a  number  of 
different  damping  and  move  limit  settings.  Instead,  they  repeatedly  violated  the  volume  constraint,  and 
the  maximum  thickness  change  diverged  as  the  iterations  progressed. 
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Figure  48:  Bending  Structure  for  Birdwing  at  Point  2 


The  maximum  thickness  change  between  iterations  of  the  two  solutions  that  were  obtained  oscillated  in 
very  high  increments  and  eventually  did  drop  beneath  the  0.01  threshold.  However,  such  unsteadiness  in 
the  tracking  data  is  indicative  of  results  that  can  hardly  be  said  to  have  converged.  In  light  of  the  poor 
solution,  Figure  48  reveals  two  separate  beam  structures  at  the  leading  and  trailing  edge  of  the  wing.  The 
proportion  and  location  of  the  “stiffeners”  at  the  end  of  the  beam  extending  towards  the  tip  of  the  wing, 
when  viewed  as  a  whole,  perhaps  resemble  feather-like  protrusions. 


As  in  the  combined  loading  case  of  Point  1,  the  distributed  viscous  drag  again  leads  to  combined 
structures  identical  to  the  isolated  membrane  structures,  and  is  shown  in  Figure  49.  When  the  viscous  drag 
is  again  removed,  the  combined  loading  case  for  Point  2  results  in  Figure  50.  In  this  case,  Point  2  elicits  a 
clearly  bending-dominated  structure.  Even  though  the  lifting  load  is  significantly  diminished,  the  induced 
drag  loads  are  extremely  low  and  hence  show  virtually  no  influence  in  the  structure,  save  a  few  faint 
battens.  Instead  the  structure  is  mainly  comprised  of  relatively  equally  sized  and  spaced  branches 
extending  from  the  root  chord  to  the  tip  chord,  like  tines  on  a  fork.  The  structure  bypasses  the  trailing 
edge  entirely. 


As  the  birdwing  in  the  perching  maneuver  begins  its  ascent,  it  is  transitioning  from  bending  loads  due  to 
lift  to  bending  loads  due  to  drag.  At  this  transitional  stage,  axial  and  normal  body  forces  are  the  same 
order  of  magnitude,  and,  more  so  than  at  the  other  points,  a  hybrid  membrane-bending  solution  is 
anticipated.  The  results  of  the  compliance  minimization  for  Point  3  are  shown  in  Figure  51  through  Figure 
53. 


The  layout  and  features  of  the  membrane  solution  in  Figure  51  are  very  similar  to  that  of  Figure  43,  with 
one  major  distinction.  The  trailing  edge  is  better  supported  than  the  leading  edge.  This  could  in  part  be 
due  to  the  forward  sweep  of  the  trailing  edge  (unlike  the  relatively  symmetrical  planform  of  the  swept 
back  configuration).  However,  Figure  51(a)  is  particularly  perplexing,  as  the  leading  edge  is  very  faintly 
supported  by  a  very  long  leg  of  a  wishbone.  One  explanation  may  lie  in  the  fact  that  there  is  less  positive 
axial  force  than  that  of  the  back  swept  configuration,  and  hence  the  trailing  half  of  the  wing  is  being 
pushed  relatively  more  in  the  negative  axial  direction  compared  to  how  much  the  leading  edge  is  being 
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pushed  in  the  positive  axial  direction.  At  any  rate,  the  disparity  between  the  support  along  the  leading  and 
trailing  edges  vanishes  with  higher  volume  fractions. 


(a)  Volf  =0.2 


(b)  Volf  =0.3 


(c)  Volf  =0.4  (d)  Volf  =0.5 

Figure  49:  Combined  Structure  for  Birdwing  at  Point  2 
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(c)  Volf  =0.4  (d)  Volf  =0.5 


Figure  50:  Combined  Structure  for  Birdwing  at  Point  2  without  Viscous  Drag 
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(a)  Volf  =  0.2 


(b)  Volf  =  0.3 


(c)  Volf  =  0.4 


(d)  Volf  =  0.5 


Figure  51:  Membrane  Structure  for  Birdwing  at  Point  3 

The  bending  structures  are  very  similar  to  those  of  the  back  swept  configuration  at  Point  1 .  A  single  beam 
extends  from  the  root  chord  most  of  the  way  to  the  tip  of  the  wing.  Also,  the  same  types  of  stiffeners 
appear.  However,  the  beam  in  this  instance  is  clearly  building  up  the  leading  edge,  and  is  not  centered 
near  the  quarter  chord. 

Figure  53  does  not  reveal  a  hybridization  of  membrane  and  bending  elements  as  anticipated.  The  solution 
again  favors  a  membrane  structure,  both  with  and  with-out  the  presence  of  the  viscous  drag.  However,  in 
the  absence  of  viscous  drag  (Figure  54)  a  thicker  leading  edge  emerges,  and  mass  starts  to  cluster  around 
the  upper  left  corner  of  the  wing,  as  it  does  slightly  in  Figure  53(d). 
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Figure  52:  Bending  Structure  for  Birdwing  at  Point  3 


When  this  happens,  the  cluster  assumes  the  role  of  a  hub,  where  separate  chains  of  truss  members  revolve 
around  the  hub  and  are  connected  to  the  hub  by  “spokes”.  The  structure  essentially  forms  a  polar  grid,  and 
consequently  little  material  is  left  to  form  battens.  The  fact  that  both  the  induced  drag  and  viscous  drag 
are  at  peak  values  accounts  for  greater  membrane  influence  and  dominance  over  bending. 

The  forward  swept  flaring  configuration  of  the  wing  is  shown  in  Figure  55  through  Figure  57  for  the  final 
point  in  the  perching  maneuver.  As  mentioned  in  the  previous  section,  the  intent  of  Point  4  was  to  capture 
the  high  braking  drag  loads  that  ultimately  deflect  the  wings  out-of-plane,  but  the  point  was  selected  when 
only  twenty  percent  of  the  velocity  remained.  Therefore  the  drag  that  should  have  produced  a  highly 
bending-dominated  structure  was  lacking,  and  instead,  the  selected  Point  4  is  actually  the  only  point 
where  the  axial  forces  (in  this  case,  due  to  lift)  exceed  the  normal  body  forces. 
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(a)  Volf  =0.2 


(b)  Volf  =0.3 


(c)  Volf  =0.4  (d)  Volf  =0.5 

Figure  53:  Combined  Structure  for  Birdwing  at  Point  3 


(c)  Volf  =0.4  (d)  Volf  =0.5 

Figure  54:  Combined  Structure  for  Birdwing  at  Point  3  without  Viscous  Drag 
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(a)  Volf  =  0.2 


(c)  Volf  =  0.4 


(d)  Volf  =  0.5 


(b)  Volf  =  0.3 


Figure  55:  Membrane  Structure  for  Birdwing  at  Point  4 

In  the  membrane  solution  at  Point  4,  the  fundamental  wishbone  structure,  while  not  completely  absent, 
does  not  best  describe  the  layout.  Rather,  three  main  fingers  (for  volume  fractions  of  0.4  and  0.5)  and  one 
thinner  finger  emanating  from  the  trailing  edge  support  the  leading  edge.  The  three  fingers  have  the 
appearance  of  three  arms  branching  out  on  a  candelabra.  They  could  alternatively  be  likened  to  the  three 
fingers  found  in  bat  wings.  Many  battens  are  again  present.  Additionally,  some  intermediate  thicknesses 
are  hanged  like  tinsel,  covering  the  bottom  half  of  the  wing,  in  the  twenty  and  thirty  percent  volume 
fractions. 

The  bending  topologies  do  not  present  any  new  features  and  are  very  similar  to  those  occurring  at  Point  1 
and  Point  3. 

In  the  absence  of  significant  drag  loading,  the  combined  structure  is  completely  membrane-dominated. 
The  combined  structure  is  again  virtually  indistinguishable  from  the  membrane  case,  and  thus  the 
combined  structure  with  no  viscous  drag  is  shown  in  Figure  57.  In  the  absence  of  viscous  drag,  a  slightly 
different  membrane  structure  emerges.  Though  there  is  no  “hub”  resulting  from  a  beam,  the  members  are 
again  arranged  in  the  fashion  of  a  web.  This  is  fundamentally  distinct  from  a  typical  truss  structure  which 
uses  the  inherent  stiffness  in  triangulation.  The  solution  is  very  clean  and  virtually  no  battens  or  stiffeners 
are  present. 
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Figure  56:  Bending  Structure  for  Birdwing  at  Point  4 


In  this  section,  both  the  structural  layout  due  to  membrane  loads  and  bending  loads  were  studied 
independently,  such  that  in  the  actual  combined  loading,  one  could  identify  whether  the  structure  assumed 
an  arrangement  in  the  membrane  or  bending  fashion.  The  bending  loads  can  produce  significantly  more 
stress  in  a  structure  than  membrane  loads  of  the  same  magnitude.  Yet,  the  combined  loading  optimization 
continually  yielded  membrane  layouts. 
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(a)  Volf  =  0.2 


(c)  Volf  =  0.4 


(b)  Volf  =  0.3 


(d)  Volf  =  0.5 


Figure  57:  Combined  Structure  for  Birdwing  at  Point  4  without  Viscous  Drag 


When  the  viscous  drag  was  removed,  the  effect  was  a  structure  that  was  a  hybrid  between  the  truss-like 
membrane  structures  and  the  beam-like  bending  structures.  Thus,  the  viscous  drag  forces  drove  the 
solutions.  Of  course,  the  normal  loads  are  still  supported,  albeit  they  were  greater  in  magnitude  than  the 
axial  loads.  So  the  implication  is  that  the  membrane  structure  with  discrete  elements  is  also  adequate  to 
support  out-of-plane  loads.  Thus  the  discrete  members  function  like  spars  in  resisting  bending,  but,  rather 
than  being  arranged  parallel  to  the  leading  edge  and  evenly  spaced,  they  are  aligned  in  a  manner  that  also 
resists  in-plane  stretching.  The  implication  of  using  a  structure  with  discrete  members  for  wing  design  is 
that  the  mechanism  designed  to  achieve  a  planform  shape  change  should  be  suitable  to  support  out-of- 
plane  loads  provided  that  the  members  have  bending  stiffness.  In  this  case,  the  actual  layout  of  the 
members  would  be  determined  by  a  mechanism  design  and  not  by  compliance  minimization.  However,  if 
compliance  minimization  is  used  to  suggest  structural  configurations,  the  general  layout  of  Points  1  and  3 
(which  are  very  similar  to  each  other)  would  probably  be  the  best  design  point,  since  the  MAV  would 
spend  the  majority  of  its  flight  time  in  those  configurations,  and  because  these  points  experience  the 
highest  loads.  Whether  a  shape-changing  mechanism  or  a  compliance  minimization  is  used  to  design 
wing  structure,  the  topologies  explored  here  demonstrate  that  the  leading  edge  must  be  predominately 
supported.  Having  a  larger  member  extend  across  the  leading  edge  may  integrate  well  with  a  mechanism 
design  where  the  leading  edge  is  simply  swept  during  a  planform  area  change. 
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3.3  Grid  Independence 

In  Section  1.2,  it  was  mentioned  that  the  discretization  scheme  for  a  distributed  parameter  system  will 
affect  the  solution,  and  thus  grid  quality  may  be  a  concern.  To  briefly  address  this  concern,  a  quick  grid- 
independence  study  is  undertaken  here.  In  a  first  attempt,  the  mesh  of  the  0°  sweep  wing  configuration 
was  refined  by  a  factor  of  1.5  and  then  1.25,  but  both  refinements  failed  in  Tornado.  It  may  be  that  the 
corresponding  cells  were  physically  too  small  for  successful  computation  by  the  vortex  lattice  method. 
However,  a  coarse  40  x  20  mesh  was  created  and  the  compliance  optimization  was  performed  for  the 
combined  loading  case  without  viscous  drag,  shown  in  Figure  54  for  the  60  x  30  mesh.  The  results  are 
shown  in  Figure  58. 

For  the  twenty  percent  volume  fraction,  little  to  no  difference  is  observed.  The  thirty  percent  volume 
fraction  case  produced  the  same  basic  layout,  but  the  coarse  mesh  lacked  one  of  the  “spoke-like” 
members  present  in  the  finer  mesh.  In  the  forty  percent  volume  fraction  case,  the  entire  outer  radius  of 
truss  members  is  missing,  and  is  not  present  until  the  fifty  percent  volume  fraction  case.  Though  the 
overall  layouts  are  of  the  same  essence,  a  finer  mesh  allows  material  to  divide  into  smaller  members, 
which  is  a  similar  effect  as  simply  increasing  the  volume  fraction.  Convergence  data  is  also  provided  in 
Figure  59. 


Figure  58:  Combined  Structure  for  Birdwing  at  Point  3  with  Coarse  Mesh 
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Figure  59:  Convergence  Data  for  Figure  58 


3 .4  Optimized  Actuation  System 

Figure  60  presents  results  of  a  wing  that  morphs  into  two  different  configurations  while  minimizing  the 
amount  of  actuation  required.  It  is  important  to  note  that  the  membrane  element  was  not  yet  implemented. 
This  explains  the  perceived  disconnect  in  some  of  the  load  paths.  Notice  how  the  optimizer  develops  the 
main  spar  along  the  leading  edge  to  support  most  of  the  load  and  then  places  various  re  volute  joints  and 
telescoping  members  to  provide  the  shape  change  actuation.  In  Figure  60  the  red  members  denote  a  “pull” 
actuator  whereas  the  blue  members  denote  a  “push  actuator.” 


Figure  60:  Results  without  Membrane  Elements 

Figure  61  presents  results  with  the  membrane  element  added.  Notice  the  increase  in  structure  to  support 
the  internal  resistance  of  the  pre-strained  membrane.  Also  the  optimizer  designs  scissoring  mechanisms  to 
achieve  the  desired  morphed  shape. 
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Figure  61:  Results  with  Membrane  Elements 
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Section  4.  CONCLUSIONS 


4.1  Summary 

In  2003,  DARPA  began  the  MAS  program  with  the  aim  of  developing  morphing  aircraft  structures  that 
can  substantially  change  shape  to  adapt  to  changing  mission  environments.  Morphing  wings  with 
planform  capabilities  of  telescoping,  chord  extension,  and  variable  sweep  will  enable  a  multi-role 
platform  far  superior  to  both  that  of  conventional  fixed-wing  aircraft,  and  even  to  aircraft  with  the 
capacity  to  modify  their  airfoil  shape.  Achieving  radical  shape  change  will  require  design  that  integrates 
innovative  technology,  including  advanced  flexible  skin  materials,  efficient  actuators  and  power  systems, 
and  mechanisms  and  structures  with  the  facility  to  execute  shape  change. 

Because  morphing  wing  design  is  a  multidisciplinary  challenge,  new  design  methods  and  processes  must 
be  synthesized  that  consider  the  holistic  problem.  The  developing  field  of  topology  optimization  promises 
to  be  an  effective  tool  that  can  simultaneously  optimize  the  arrangement  and  connectivity  of  the  load- 
bearing  structure,  target  shape-matching  mechanism,  and  flexible  skin  elements,  for  any  number  of  design 
objectives.  Therefore,  two  objectives  were  established  to  advance  previous  efforts.  The  first  objective  was 
to  develop  a  six  degree-of-freedom  finite  element  model  comprised  of  continuous  elements.  This  model 
would  then  be  capable  of  supporting  out-of-plane  loading  and  could  be  used  to  represent  the  flexible  skin 
of  the  wing.  It  was  important  to  implement  the  model  in  MATLAB  to  be  able  to  easily  integrate  it  with 
optimization  routines.  The  second  objective  was  to  consider  realistic  aerodynamic  loading  that  MAV 
wings  might  experience  through  a  changing  mission  environment.  The  perching  maneuver  common  to 
birds  was  selected  as  an  interesting  miniature  of  a  multi-role  mission  in  which  a  bird  undergoes  several 
planform  transformations. 

By  constructing  a  compliance  minimization  optimization  objective  (equivalent  to  a  stiffness  maximization 
objective),  conceptual  wing  structure  layouts  of  a  MAV  at  various  points  throughout  a  perching  maneuver 
were  investigated.  Both  the  structure  corresponding  to  the  isolated  cases  of  membrane  and  bending 
loading,  as  well  as  the  combined  loading,  was  optimized.  The  results  revealed  some  distinctive  features  of 
wing  structure  when  only  membrane  loads  or  when  only  bending  loads  were  considered.  The  structures 
formed  by  membrane  loads  were  comprised  of  many  discrete  truss-like  members  that  generally  supported 
the  leading  edge  of  the  wing.  Rather  than  predominantly  relying  on  triangulation  for  inherent  strength,  the 
members  were  typically  arranged  in  cobweb  formations  with  quadrilateral  voids.  Straight  rods  or  battens 
protruded  from  the  truss  members  pointing  towards  the  center  of  a  voided  region  or  towards  the  outskirts 
of  the  wing.  Contrary  to  the  membrane  structures,  the  bending  loads  elicited  a  principal  beam-like 
structure  situated  anywhere  between  the  leading  edge  and  the  quarter-chord,  and  extending  from  the  root 
chord  to  the  tip  of  the  wing.  With  additional  material,  the  shape  of  the  beam  evolved  into  a  wishbone  or 
delta  conglomeration.  Contrary  to  the  battens  that  formed  in  the  membrane  case,  chutes  with  endings  that 
branched  in  arbitrary  directions  emanated  from  the  central  beam.  These  chutes  act  like  stiffeners  similar  to 
the  neural  structures  found  in  the  flat  wings  of  insects.  The  full  three-dimensional  loading  typically 
favored  the  membrane  structural  formation  over  that  of  the  bending.  However,  in  the  absence  of  the 
viscous  drag  distributed  over  the  surface  of  the  wing,  a  structure  formed  that  is  a  hybrid  of  the  typical 
membrane  and  bending  structures.  Thus  the  truss-like  structure  corresponding  to  the  combined  loading  is 
thought  to  be  driven  by  the  distributed  viscous  forces.  However,  the  discrete  members  must  also  behave 
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like  spars  that  can  resist  the  bending  loads,  which  ultimately  contribute  the  most  to  the  total  magnitude  of 
the  forces  on  the  wing. 

4.2  Recommendations 

Though  the  resulting  structures  studied  in  Section  3  revealed  distinctive  features  that  parallel  many 
formations  with  inherent  strength  found  in  nature,  the  methodology  developed  in  Section  2  can  certainly 
undergo  improvements.  First  of  all,  the  discretization  scheme  of  the  geometry  for  the  finite  element  model 
was  chosen  to  be  structured  meshing.  Though  the  structured  mesh  was  much  easier  to  implement  and  also 
allowed  a  common  mesh  to  be  used  for  both  the  finite  element  method  and  the  aerodynamic  vortex  lattice 
method,  an  unstructured  mesh  is  more  versatile  in  accurately  modeling  complex  geometries.  The  inability 
of  the  structured  mesh  to  model  highly  non-rectangular  geometries  was  observed  with  the  dive 
configuration,  and  merits  greater  consideration  of  using  an  unstructured  mesh.  However,  if  an  un¬ 
structured  mesh  is  to  be  employed,  the  advantage  of  easily  discretizing  a  geometry  with  triangles  is  lost, 
since  the  element  geometries  derived  were  quadrilaterals. 

The  finite  element  itself  could,  of  course,  be  selected  differently.  For  example,  rather  than  using  Kirchhoff 
plate  theory,  Mindlin-Reissner  plate  theory  could  be  used  to  couple  the  in-plane  and  out-of-plane  effects. 
This  would  result  in  a  five  degree-of-freedom  element  that  would  avoid  the  need  for  superimposing  two 
separate  membranes  and  bending  elements  together.  Such  an  element  should  also  avoid  the  small 
displacement  assumption  of  the  combined  membrane -bending  element.  Also,  higher-order,  nonlinear 
shape  functions  could  be  used,  rather  than  the  simple  bilinear  shape  functions  used  here.  Higher-order 
shape  functions  should  increase  the  accuracy  of  the  finite  element  model,  particularly  for  cambered 
geometry. 

The  MATLAB  implementations  were  very  simple  and  the  routines  ran  very  efficiently.  Both  the 
membrane  and  bending  models  were  validated  with  classical  solutions,  and  their  accuracy  was  verified. 
The  superimposed  membrane  and  bending  element,  while  not  validated  with  an  analytical  solution,  did 
match  the  individual  solutions  of  the  membrane  and  bending  models.  The  three  dimensional  six  degree- 
of-freedom  model,  however,  still  lacks  validation  and  successful  integration  with  the  SIMP  model  of  the 
compliance  optimization.  If  a  full  three  dimensional  analysis  is  not  easily  obtainable  for  comparison  and 
validation  of  the  model,  the  three  dimensional  FE  model  can  be  validated  with  a  commercial  finite 
element  code,  such  as  NASTRAN. 

The  aerodynamic  force  estimations  can  afford  much  improvement,  though  the  extent  to  which 
improvement  would  influence  the  optimized  thickness  distributions  may  not  be  substantial.  Effort  should 
be  spent  to  ensure  that  the  outputted  forces  from  Tornado  do  in  fact  act  through  the  centroid  of  each 
panel.  Post-stall  aerodynamic  models  can  also  be  employed  to  improve  the  aerodynamic  calculations  at 
the  end  of  the  maneuver.  Also,  a  more  ideal  selection  of  Point  4  in  the  perching  maneuver  would  capture 
the  high  angle  of  attack,  and  high  angle-of-attack  and  pitch  rates,  while  the  aircraft  still  has  most  of  its 
velocity.  This  should  in  turn  capture  the  high  bending  loads  that  the  drag  exerts  on  the  aircraft  as  it  begins 
to  slow  down  into  a  vertical  landing. 

In  general,  the  Optimality  Criteria  method  performed  very  efficiently,  usually  completing  a  minimization 
in  less  than  five  minutes  with  a  couple  seconds  at  most  between  iterations.  The  Method  of  Moving 
Asymptotes  (MMA)  is  often  employed  for  large  scale  optimization  problems  and  may  reduce 
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computation  time.  If  larger  problems  are  undertaken,  such  as  would  be  produced  from  mesh  refinement, 
the  MMA  solver  may  be  worthy  of  investigation.  However,  for  the  problem  sizes  explored  in  Section  3, 
the  OC  method  executed  very  well. 

4.3  Future  Work 

With  the  insight  gleaned  from  structural  layouts,  some  experimentation  is  probably  in  order  to  test  the 
stiffness  of  the  structures.  For  instance,  a  study  could  be  done  to  determine  how  much  rigidity  battens  add 
to  a  structure  in  two-dimensional  stress,  or  how  much  rigidity  veins  or  branches  add  to  a  beam 
undergoing  bending  stress.  Also,  the  inherent  strength  of  a  wishbone  structure  to  resist  being  spread  apart 
could  be  explored.  A  final  pattern  observed  in  the  results  had  a  layout  similar  to  a  cobweb.  To  aid  the 
investigation,  these  structural  formations  could  be  more  directly  correlated  to  the  wings  of  birds  and 
insects. 
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LIST  OF  ACRONYMS,  ABBRECIATIONS,  AND  SYMBOLS 


AFRL 

Air  Force  Research  Laboratory 

AIC 

Aerodynamic  Influence  Coefficient 

DARPA 

Defense  Advanced  Research  Projects  Agency 

DOF 

Degree  of  Freedom 

FE 

Finite  Element 

GUI 

Graphical  User  Interface 

KKT 

Karush-Kuhn  Tucker 

MAS 

Morphing  Aircraft  Structure 

MAV 

Micro  Air  Vehicle 

MMA 

Method  of  Moving  Asymptotes 

OC 

Optimality  Criteria 

SIMP 

Solid  Isotropic  Material  with  Penalization 

SMP 

Shape  Memory  Polymer 

UAV 

Unmanned  Aerial  Vehicle 

VLM 

Vortex  Lattice  Method 
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